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SUMMARY 


The steady axially symmetric stagnation point flow of an incompressible electric- 
ally conducting viscous fluid in the presence of a magnetic field normal to the wall 
is investigated. The wall is assumed to be thermally insulated and all physical 
properties of the fluid such as viscosity, electrical conductivity, and magnetic 
permeability, etc., are assumed to be independent of the temperature and the strength 
of the magnetic field. Solution of the equations of magnetohydrodynamics leads 
to the conclusion of the existence of three regions of flow: (i) a layer of fluid near the 
wall in which viscosity is important, called the magneto-viscous layer having thick- 
ness of order (v/a)', (ii) a buffer-layer in which viscosity is unimportant, called the 
magneto-inviscid layer and having thickness of order 1/(87op, a)*, and (iii) a region 
of potential flow in which the velocity field is directly proportional in magnitude to 
the magnetic field. Here a is a parameter of dimensions (time)~!, such that the 
radial component of external velocity is ar. 

A numerical solution of the basic equations has been obtained for the special 
ease of 8B — A2v = 10° and P = 0-72, in the form of a series expansion in terms 
of the magnetic parameter m pe? Hi(2c)/(pa). In particular it appears that the 
presence of the magnetic field produces, in the vicinity of the stagnation point, 
a considerable reduction in the local shear stress and the eigentemperature at the 
wall. The new results obtained are intended to supplement the results of Neuringer 
and McIlroy (1), Rossow (2), and Meyer (3) who have investigated the correspond- 
ing problem of the two-dimensional stagnation point flow. 


Symbols 


E.M.U. and C.G.S. system used throughout. 
electric intensity. 
(H,,0,H,) magnetic intensity. 
current density. 


electrical conductivity. 
magnetic permeability. 
(u,,0,u,) fluid velocity, 
fluid density. 
dynamic viscosity. 
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p/p kinematic viscosity. 
parameter of dimensions (time)~'. 
pressure. 

temperature. 

thermometric conductivity. 
specific heat at constant pressure. 
viscous heat dissipation. 

,© dimensionless functions describing the local temperature distri- 
butions in the magneto-viscous and magneto-inviscid layers 
respectively. 

Prandtl number. 
| 4rou, magnetic viscosity. 
| S7ou,v a parameter. 
pw? H2 (2c) pa magnetic parameter. 
gg, h dimensionless functions of € = (a/v)!z relevant to the magneto- 


viscous layer. 
 (, H dimensionless functions of = (8zop,a)'(z-+-A) relevant to 


the magneto-inviscid layer. 
(r,d,z) evylindrical polar coordinates. 


Subscripts 
w body surface. 


~ free stream. 


1. Introduction 
RECENT papers by Neuringer and Mellroy (1), Rossow (2), and Meyer (3) 
deal with the reduction of shear stress and heat transfer rates due to the 
influence of a magnetic field on two-dimensional stagnation point flow. 
The present paper deals with the magnetohydrodynamic effects on axially 
symmetric stagnation point flow and heat transfer. For simplicity, as in 
references (1) to (3), we shall consider a viscous incompressible electrically 
conducting fluid, and all physical properties of the fluid such as viscosity, 
electrical conductivity, magnetic permeability, ete., are assumed to be 
independent of the temperature and the strength of the magnetic field. 
Cylindrical polar coordinates (r, ¢, z) are used, the wall coincides with the 
plane z = 0, and the stagnation point is at the origin. Also 

(a) the longitudinal flow is taken to be in the negative z-direction, 

(b) the magnetic field at the wall is taken normal to the wall, and 

(c) the wall is thermally insulated, i.e. (é7'/éz),_9 = 0. 
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The basic equations for axially symmetric flow in E.M.U. and C.G.S. 
units are, in the usual notation: 
divq = 0, (1.1) 
Dq 1 . , 
9.” grad p = vV2q+"¢jx H, (Lt 
Dt p p 


: a es “wir, (1. 
together with the Maxwell electromagnetic field equations, 
divH = 0, 
curlH = 47j, 
cH 
em’ 


curlE = —p 


and finally Ohm’s law for the moving fluid 
j — o(E+p,qH). (1.7) 
In equation (1.5) the displacement current is neglected; for steady motion, 
however, the equation is exact. 
For steady motion, on taking the curl of equation (1.7) and using (1.6) 
and (1.5), there results 
V x(VxH) = 47op, V x (q * H). (1.8) 
On expanding V x (V x H) and V x (q x H) and using equations (1.1) and 
(1.4), equation (1.8) becomes 
AV*H = (q.V)H—(H.V)q (1.9) 
where A — 1 47op,; A is defined as the magnetic viscosity. Also on using 
equation (1.5) the equation of motion for steady flow becomes 


TT, 


lV p4 (q.V)q = vV2q +e. (Vx H) x H. (1.10) 
p 4p 


In cylindrical polar coordinates equations (1.9), (1.10), in component 
form, and equations (1.1), (1.4), give the set of equations: 


= (ru,) + <(ru,) = 0, (1.11) 
. Ps, lé pour , eu, ] ben eH, _¢H, ’ 
p or rér\ or} é2% rr] dap “\ez or 
(1.12) 
; ae < le yous , Ou, 4 Me eH, _H, 
p oz rér\ er] @2%| 4p "\or  & 


(1.13) 
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“ (rH,) + (rH) = 0, (1.14) 


CT c 
; L 2(-2t % 
ror cr 


yf! (rth) " fey (1.16) 


cr cr ror cr 
Finally, the energy equation (1.3), which describes the transport of 
thermal energy in the presence of a magnetic field is, for axially sym 
metric steady motion, 
oT of Leal et @ - ; 
pe,(u, 1 2, -k ( }. | ++. 3% (1.17) 
cr ~ OB ror\ -cor 


where ® is the viscous dissipation function, and has the value 


Gu,\?  feu,\? u? feu,  eu,\? 
(p Zu | " +( ‘| } r } coe - 
cr oz go Bh oz cr 
The term (1/c)j? is the Joule dissipation, and 


7? i (1.19) 


aoe 
oO oO l62*0 


ho. l (“t =) 


cr C2 


2. Similarity solution of the equations, and boundary conditions 

Consider now the magnetohydrodynamic equations (1.11) to (1.16). 
The boundary conditions appropriate to the flow under consideration are 

u u. 0 for z 0 
iii l (2.1) 
and u,-ar as 2>@ forallr 
where a is constant. Following Homann (4), who discussed the similar 
solution for the field free case, we introduce the new variable 
(a/v)*z, 

and the velocity components 


q (u,, O, uw.) |raf’(f), 0, — 2(av)*f(C)}. (2.3) 

When these substitutions are made in the basic equations of motion 

(1.11) to (1.16), a similar solution for axially symmetric stagnation point 

flow in the presence of a magnetic field is found to be possible if the law 
of deformation of // in component form is taken as 


a\} 
H = (//,,0, 7.) (—1,(“)*ri'(Z),0,241,.4(0)) (2.4) 
> 


where //,, is a representative magnetic intensity. The above expressions 
for q and H satisfy identically equations (1.11) and (1.14), and equations 
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(1.12) and (1.13) become: 


ra*( f”" +-2ff" —(f')? —2ahh")— . 
p 


0, 
cr 


* 
va( aff’ +-2f")+ aa*rh’h" + 5 


where « — pu, H?./42pva. The pressure may be eliminated between equa- 


0, 


tions (2.5) and (2.6); the resulting equation can be integrated once with 
respect to ¢, and we obtain an ordinary differential equation defining f as 
follows: | i 2ff’"—(f')?? n x{ (h’)?- : 2hh" | | C 0, (2.7) 
where C is a constant of integration. The boundary conditions on f are 

f(0) = f'(0) = 9, f'(o) = 1. (2.8) 
The corresponding equations for the magnetic field, namely (1.15) and 

6), are ” ” ” 

(1.16), are now Bh hf’ —fh’, (2.9) 
and Bh” = hf'—fh’, 2.10) 
where £ is the ratio of the magnetic viscosity to twice the kinematic 
viscosity. Since equation (2.10) is a first integral of equation (2.9) the 
validity of the similar solution assumed in expressions (2.3) and (2.4) is 
established. It remains now to state the boundary conditions for H. The 
similar solution for H implies, from the usual condition on B and H at 


the solid interface, that the field in the solid must satisfy the boundary 
conditions 


? : 
H — (H,, Mp, H,) (—1.(")'rH(, 0, 2yH, (0) 2.11) 
: 


at z — 0 for all r, where y is the ratio of the magnetic permeability of 
the fluid to that of the solid. If we suppose the solid is non-conducting 
the magnetic field inside the solid is determined by H = Vd, where ¢ 
satisfies the Laplace equation V*¢ = 0. If now we impose the restriction 
that the magnetic field is to remain finite as z-» —® then in the solid 
the only possible solution of the form required is 


H, = Hy= 9, and H, = 2yH,,h(0). 2.12) 

Thus for the fluid at the solid interface the boundary conditions on H are 
H (H,, Ho, H.) (0,0, 2H) and hence 

A(O) 1, and h’(0) = 0. (2.13) 

For a physical interpretation of the above situation in which the mag- 

netic field in the solid is uniform and in the z-direction we might consider 

the whole assembly to be placed inside a large solenoid which provides 

an applied magnetic field in the solid, so that the field is approximately 
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uniform in the neighbourhood of the axis r =~ 0, and acts in the z-diree 
tion; the resultant field in the fluid is the sum of this field and that due 
to the motion of the fluid. A model corresponding approximately to the 


solution discussed here is that in which the whole assembly becomes 
infinitely large.4 


Returning to equation (2.7) the undetermined constant C may be 


specified by examination of the solution of equations (2.7) and (2.10) when 
Cis large, Since, for large (, f'(0) ~ Land f(t) ~ ¢+ A, then from equation 
(2.10) A°(0) ~ Band hil) ~ Bit +A), where A and Bare further constants 
of integration. Furthermore A°(0) ~ © and if we seek a solution for f in 


which f"(¢) and all higher derivatives vanish identically, then 
C = 1—a* = 1—afh'(0)))_.. (2.14) 


Hence, to recapitulate, we need to solve the following coupled system of 
ordinary differential equations 


2ff" 4+-1—(f')* (h')?— 2hh"—(h')y? | 2.15) 
i -+-fh' —f ‘he ), (2.16) 

subject to the boundary conditions 
iT 0, f’lco) . Alo) . ante) 0. (2.17) 


t, H* Anpyva and B A 21 (Srop,v) '. We note that when 
i O, equations (2.15) and (2.16) uncouple and equation 


15) reduces to the field free case discussed by Homann (4) and Fréssling 


3. Series solution of the ordinary differential equations (2.15) and 

(2.16) 

The numerical solution of equations (2.15) and (2.16) subject to the 
boundary conditions (2.17) is complicated by non-linearity, coupling, and 
the fact that for most practical purposes £ is very large, say of the order 
1 or 10%) also equation (2.15) is of the boundary value type and (2.16) 
is of the initial value type.- The complication of non-linearity and coupling 
is reduced by expressing f and A as a series in powers of the magnetic 


parameter m 1h pr 12.20 pa, tor fixed Bp; i.e, 
f(C) S mf tl), h(C) YS mh (0). (3.1) 
— r ~_ r 


r 0 r 0 


On substituting the series (3.1) into equations (2,15) and (2.16) and 


equating coeflicients of like powers of m, there results a set of differential 


* The authors are indebted to a referee for suggesting this model 





ON AXIALLY SYMMETRIC STAGNATION POINT FLOW 


equations, of which the first three are: 
S04 Yofo—(Sot+1 ~ 0, (3.2) 
Bho+foho—Soho = (3.3) 
fy + Beli —AoSs+ Aol, = Pl[2heho—(ho)* a (3.4) 

subject to the boundary conditions 

f(9) — S509) 0, Ll) : (3.2 a) 
ho) l, hi (0) 0. (3.3a) 
fo fA) 0, Si) , (3.4a) 

Consider first equation (3.2). From the numerical solution for fo, 
obtained by Fréssling (5), we observe that for ( -= (y (where G, — 3:7) 
we have f, —§ ©» A. Thus there are two regions of interest: (a4) an “inner 
region’ for © - G in which viscosity is important, where f, and its higher 
derivatives vary rapidly from their initial values at (0 to their ‘asymp 
totic values’ at © %, and (4) an ‘outer region’ of flow in which the 
Viscosity is unimportant, and where f, © (4A for (<< 0 << ®. 

Equation (3.3) and the boundary conditions (3.3a) define an initial 
value problem for hy. It is convenient to transform equation (3.3) by 
introducing the new dependent variable /,, defined by 

h , 
hy 14 ‘0 (3.5) 
B 


so that the equation and boundary conditions become 


hi, fot (hols Kifod, = gl) — ho) — 9. (3.6) 


The numerical solution of this equation may be obtained by step-by-step 


methods.t When ¢ > {, equation (3.3) becomes 


where for clarity the outer solution for A, is denoted by //,. Hence 


HI, hy + col’ 47" 7] | e dy) (3.8) 
” 

where .By — C+ A, and by and ¢, are constants of integration which may 
be found by joining the inner and outer solutions at ( — (). The ‘inner’ 
and ‘outer’ solutions to hy defined by ha(Q) and H,() respectively are now 
known to the same accuracy as the accepted tables of f,, and in this sense 

we can say hy is ‘exact’. 
+ Note that when 8 10° there is fortuitous cancellation inside the bracket on the right- 


_ 
hand mde of equation (3.6). Hence a good approximation to hy is {fog for 0 << <= %. 
0 
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The equation for f, is of the ordinary boundary value type and we can 
solve this equation using step-by-step methods, bearing in mind that we 
must consider the ‘inner’ and ‘outer’ solutions of f, corresponding to the 
‘inner’ and ‘outer’ solutions of hy. If the operator L(D) is defined by 


3 ° d? or d 1 ofr 
L{D) (an 2fo ie Fe wet fo): 


then the ‘inner’ solution for f, is obtained by solving 


L(D)f, = 624 2(1 ; Ko 


For ¢ < (, the required solution of equation (3.9) may be written 
f=SfPt+nJsi, (3.10) 


where f\”’ is a particular integral and f\? is a complementary function, 
each of which may be chosen to satisfy the boundary conditions on f, at 
¢- 0, The constant y, remains to be determined from the boundary 
conditions at infinity. 


When ¢ Cy, equation (3.4) becomes 
Si +204 A)fi—2f, = Bho)? . +B{2hahg—(ho)?}. (3.11) 


Since a numerical solution of this equation is necessary, it is convenient 
to make the change of variables 


(+4 = By, fl) = vBR(n), (3.12) 


so that the equation becomes 


l we ” ’ » ” , 9 . ‘ 
SFP + 2 FF) = (Hy)3 2+ 2M (Hg)? = U(n), say. (3.13) 
i 
The solution of this equation must be discussed in some detail, since the 
method applies to subsequent equations; the physical meaning of the 
approximations involved is considered later. The parameter f is very 
large, and it is reasonable to expect that a suitable approximate solution 
of the equation may be obtained by neglecting the term F7'/8, so that the 
equation becomes ” oy “a 
2(7F{— Fi) = t(n). (3.13) 
The third-order equation (3.13) is thus replaced by the second-order 
equation (3.13a), so that one boundary condition must be abandoned; 
furthermore, the solutions of both equations must be compared. In the 
first place we can obtain a particular integral F\”) which satisfies equation 
(3.13 a), and also equation (3.13) correct to terms O(1/8). 
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The general solution of equation (3.13) is then 


7 7 7 
2 ( 9-Bnt oe 1 ,-Bn'| 2a. Fe 
F, a a7 | e By dy 4 2B | e By dn + 59° Bn | +A, +d3+ Fy ) 


. 


0 0 


(3.14) 
and that of equation (3.13 a) 
, 2 i ¥( ) J 
F, = 6, ?+6,+ FY’. (3.14a) 
These functions differ by the term in a,, and since we require to join the 
inner and outer solutions for f, at ¢ — C = 3-7 we consider the behaviour 
of this term for » > 7, where mg is the corresponding value of 7. 
Asymptotic expansions may then be used for the integral, and the coeffi- 
cient of a, is found to be 


La ihe yerrtae 


For the range of values of » under consideration, the first two terms alone 
in the above expression are significant (the exponential terms are of order 
e~*), so that the solution (3.14) is, effectively, the same as (3.14a). Since 
the particular integral F\” is correct to order 1/8, we may expect this 
accuracy in the evaluation of F,. 

The boundary condition at infinity is satisfied by choosing 6, so that 


d 
In 


and the value of y, in (3.10) is then obtained by ensuring continuity in 
df, df at € — %; by is obtained by numerical integration. The results 
show that the term neglected in equation (3.13) is indeed small. 

The higher order functions h,, f,, 4,, and f,; may now be obtained in a 
similar fashion, i.e. we calculate the ‘inner’ solutions defined by f,(¢) and 


B 


where 4° is the Kronecker delta, together with the ‘outer’ solutions defined 
by F(») and H,(y) respectively. We note that the ‘inner’ and ‘outer’ 
solutions for f,() (r > 1) are thus known approximately to within terms 
of order 1 8, since in general we have neglected the term (1/8)F,” in the 
differential equations which describe the outer solution. As this term 
originates from the viscous terms in the original flow equations (1.12) and 
(1.13) the above approximation is equivalent to stating that the outer 
solution applies to a region in which the viscosity is unimportant. In the 
next section a physical interpretation of this approximation is discussed. 


2b, n+ Fv >0 asn>D, 
( 


hl) = b+ 2h), (3.15) 
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4. The magneto-viscous and magneto-inviscid layers 

The analysis of the previous section reveals the existence of three régimes 
of flow. In the first place there is a layer of fluid near the wall, which may 
be described as the magneto-viscous layer, and this is of the conventional 
boundary layer type. That is, inertia terms and viscous forces are com- 
parable in magnitude, and the thickness is of the order (va)! as in the 
field free case. However, the Lorentz force is important in this layer, and 
furthermore the pressure distribution throughout the layer and the velocity 
distribution near the edge of the layer are greatly influenced by the 
adjacent magneto-inviscid layer. 

The magneto-inviscid layer plays the role of a buffer layer between the 
magneto-viscous layer and the potential flow régime in the free stream. 
In this layer viscosity is unimportant, but the Lorentz force is comparable 
in magnitude with the inertia terms. From section 3 we find that this 
layer has thickness of order (Sz7ou,a)~', so that if B = O(10%) the ratio 
of this thickness to that of the magneto-viscous layer is O(10%). Again, 
from section 3, the error introduced by considering this layer to be 
inviscid is seen to be of the order 8-', and is thus extremely small in the 
case under consideration. Thus equations describing the field and the 
flow in this layer may be obtained by putting v 0 in equations (1.11) 
to (1.16). If, further, we introduce for this layer the variable » defined 
in (3.12), and velocity and field components defined by 

q = |raF'(y), 0, —2a(8zop,a)-* F(m)], (4.1) 
H | —H. (Szopu,4a)'rH'(y), 0, 2H, H(n)|, (4.2) 
we obtain ordinary differential equations for F and H by the process 
described in section 2 for f and 4, Expansion in series as in section 3, 


namely, 


F(n) = S m'F(n), — H(y) = 3 mH), (4.3) 
r ) 0 


( r 
then yields a set of equations for F. and H,. 

Thus the solution described in section 3 is obtained by joining the 
solution in the magneto-viscous layer to that in the magneto-inviscid 
layer. Since the continuity of velocity and field components u,, u,, H,, H, 
must be ensured at the interface, the following conditions are obtained 


at ny = (C9+-A)/B: 


TS A(Ge) \BF ( No)s 8° + hale) - H,( yp) | 


" (4.4) 
(ie). 


Above the magneto-inviscid layer lies the free stream; this régime of 


dF. (¢ h, jaft, , 
(i). (a) =o) J 


\ ‘ne 
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flow is included in the equations for the magneto-inviscid layer, since it 
is represented by the solution for large values of y, but it is convenient 
to think of it physically as a separate layer of flow as in boundary layer 


theory. 
Tables 1-4 give the functions f,.,(¢), F..,(7), 4,(¢), and H,(y) for r = 9, 
1, and 2, together with their associated derivatives, when 8 = 10°. 


5. Solution of the energy equation when the wall is thermally 
insulated 
The partial differential equation which describes the transport of 
thermal energy for axially symmetric flow ist 


pe, (uo , “.=} k l (rt) ad 2 


cr C2 | ror\ or 0z* 


>(eu,\? , 9(euz\? , 2u2 (eu, , eu,\* | foo, O@Y 
pl 2H) oft) 2m (Ate Smal] 1 (Cll PHLY? 6.1 
Or 7 | Oz - Cz Or | | l627o0\ cr 02 
where the second and third terms on the right-hand side of equation (5.1) 
represent the viscous and Joule dissipation respectively. We shall con- 
sider the case of a thermally insulated wall, i.e. 
C T 
at z = 0, -0 for all r | 
cz 


and TT, asz->x 


where 7, is the bulk temperature of the fluid in the main stream. 

Equation (5.1) may be analysed as follows. First we determine the 
conditions for which a similar solution exists as in section 2, and then 
proceed to calculate the outer solution describing the local temperature 
distribution correct to within quantities of order 8~' as in section 3. 
Thus we find on examination of equation (5.1) that for the transport of 
thermal energy it is sufficient to consider: 

(4) a magneto-viscous thermal layer in which the convection, conduc- 
tion, viscous, and Joule dissipation terms of (5.1) are important. Further- 
more we neglect the first three terms of the viscous dissipation function 
and also the radial conduction, all of which are small in comparison with 
terms of the same order as the axial conduction. 

(B) a magneto-inviscid thermal layer in which axial and radial conduc- 
tion are negligible. Since in this region vy = 0, then only convection and 
Joule dissipation contribute to the transport of thermal energy. 


+ Note that equation (5.1) has limited application since (7’'—7',)/T,, must be small. 
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A. The magneto-viscous thermal layer 
We introduce a dimensionless temperature distribution 
, c,(T' —T, 
6(C) i ) 


=) 
a*rT- 


and substitution in equation (5.1) gives the equation 
po + 2( f0’ —f'0)+-(f")?+-2mp(h")? = 0, 


where P is the Prandt! number. One boundary condition is 
A’(0) 0, (3.35 
and the remaining condition must be obtained by ensuring continuity in 
the temperature at the interface between the magneto-viscous and mag 
neto-inviscid thermal layers. 
For small m we write 


“af)= P> 


{ 


mq,(C), (5.6) 
) 
and on using the expressions (3.1), (3.15) in (5.4), and equating the coeffi 
cients of the same powers of m, we have 
Jo +2P(fogo—fode) (fo)*, 
G+ 2P(fogi—foa) = 2PULi Gof Go) — Sofi —(ho)*, 
and so on. The boundary conditions for the g, are 
gi{O) = 0, allr. (5.9) 
The remaining set of boundary conditions are introduced in the discussion 
of the magneto-inviscid thermal layer. 
B. The magneto-inviscid thermal layer 
In the magneto-inviscid thermal layer the mechanism for energy trans 
port consists of convection and Joule dissipation, i.e. 
oT 
pe, (u,—— 4 


u (5.10) 


=| l on =). 


cr ” Ga 


l62%0\ cz cr 

With the introduction of the new variable » defined in (3.12) together 
with velocity and magnetic field components, as defined by (4.1) and (4.2), 
and also the dimensionless local temperature distribution 
c,(7'—T, ) 


(7) 2p? 
a*r? 


(5.11) 


we obtain F’'®™) — Fo’ m(H")?, (5.12) 
to be solved subject to the boundary condition 


(-)/ ”) >0 as n> 2. (5.13) 
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Since for the field free solution ©(7) = 0 when the Prandtl number is of 
order unity (for gases), we may write for small m 


(7) > m’G,(y). (5.14) 
rel 


Substitution of the series (5.14) and (4.3) into equation (5.12) gives the 


equations nG,—G, (H")2, (5.15) 
7G,—G, — FL G,—F,G,—2Hg Hj, (5.16) 

and so on, where G(») > 0 as n> &. 
Proceeding as in section 3 we now calculate a particular integral and 
a complementary function for g,(¢) in the interval 0 << (< %. The 
appropriate combination is then determined from our knowledge of the 
particular integral G,(»). (The complementary function of G,(7) is un- 
acceptable, since G(x) — 0.) In effect since the temperature must be 
continuous at the interface between the magneto-viscous and magneto- 


inviscid layers, then Pg(Lo) = Gn) (5.17) 


for r i See 

In Tables 5 and 6 the functions g,(¢) and G,(y) P for r — 0, 1, 2, and 3, 
together with their first derivatives, are given in tabular form, The Prandtl 
number was taken to be 0-72. 


6. Results and discussion 

From Tables 1-6 quantitative information may now be deduced for 
How characteristics such as the local shear stress and local temperature 
at the wall; the results are relevant only to the special case 8 — 10% and 
P — 0-72. 


Consider first the local shear stress at the wall defined by 


: 1 
2) tore 
C2 }2-9 v 


If +, is the field free local shear stress at the wall, we obtain 
T I 


=" fro) fo) + mf (0) + mf (0) + mefs(0) 4 


where f5(0) 1-311,, and from Table 1, 
fi) 2-765s, S30) — 0-368, f3(9) —I1-161. 

The ratio r,. 7) is given graphically in Fig. 2 for various values of m. 
Owing to the truncated form of the calculated series (6.2) this ratio is 
kriown with certainty only for 0 < m <— 0-2. Thus for m — 0-2 there 
is a 41 per cent reduction in the local wall shear stress. However, it is 
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instructive to observe that the local shear stress appears to be zero for 
a value of m near to 0-46. When the field is as strong as this, the resulting 
type of flow pattern is not clear. Tentatively one might expect that when 


7, = 0 a region of stationary fluid exists near the wall. The possibility 
of a solution existing for 7, < 0 (i.e. m greater than the estimated value 


0-46) in which we should have a radial inflow along the wall seems to be 
unlikely. Obviously the need for further computations in the neighbour- 
hood of m 0-4 is indicated. 


fis) Fn) 
+04 1 


08} 
06+ 
04. 


02- 





4 . + — 1 4 


’ 2 3 0 1 2 s 
Fic. 1, Dimensionless radial velocity profiles, u(r, z)/u(r, 0) = f'(f), F’(y), 
for various values of m, when B 108. 


The dimensionless radial velocity profiles f’({) for the magneto-viscous 
layer and F’(y) for the magneto-inviscid layer are given graphically in 
Fig. 1 for m 0, 05, 0-10, and 0-20, Physically the influence of the 
magnetic field on shear stress and velocity profiles may be described as 
follows. As the radial velocity in the magneto-viscous layer increases 
from zero at the wall the induced current increases and consequently the 
Lorentz force ,:,(j « H) grows in strength. Since this force has a radial 
component which opposes the radial flow, the radial velocity will be 
retarded and thus reduced in magnitude as shown in Fig. 1. The radial 
velocity component will then increase slowly through the magneto-inviscid 
layer to the incident value in the layer of potential flow, as the Lorentz 
force decays slowly to zero. If in the magneto-viscous layer the magnitude 
of u, is reduced, then so will be the local shear stress at the wall. 

The effect of the magnetic field on the local temperature distribution 
at the wall, known as the eigentemperature, may be seen on examining 
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. 2. Shear stress and eigentemperature ratios versus m when B 10° 
and P 0-72. 


2 


. 3. Dimensionless local temperature profiles, ¢,(T— T,,)/a*r? = 6({), 
©(»), for various values of m, when 8 = 10° and P = 0-72. 


the ratio of @,, = ¢,(T,,—T,)/a*r* with magnetic field to @,, = 0) without 
the magnetic field. Thus 


— i | | 
* ae Joy 900) +m gu(O) +m AGa(O) + mAgs(O) +...}, (6.3) 
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where from Table 5: 


Jo(9) = 0-588,  g,(0) —0:730, g.(0) = —O0-758;, g3(0) = —0-795. 








J 


6 


Fic. 4. Representative dimensionless magnetic field profiles (see expressions 
(2.4) and (4.2)), when B 10° and 0 < m — 0-2. 


The ratio 4, 4) is given graphically in Fig. 2 and representative dimension- 
less thermal profiles (¢,©) are given in Fig. 3 for m = 0, 0-05, 0-10, and 
0-20. We note that when m = 0-2 there is a 31 per cent decrease in eigen- 
temperature. Furthermore, assuming that equation (6.3) is valid for 
m > 0-2 we observe from Fig. 2 that the eigentemperature will be zero 
for m slightly in excess of 0-46. The question of whether or not such a 
situation is possible remains unanswered. However, the reason for the 
reduction in eigentemperature 4, and the local temperature distribution 
throughout the magneto-viscous and inviscid thermal layers can be seen 
on examination of the energy equation (5.1). In the magneto-viscous 
thermal layer near the wall the Joule dissipation will be negligible but 
there is a considerable reduction in the viscous dissipation, i.e. the term 
p(eu, €z)*. Since the energy equation (5.1) is linear in 7’ it follows that the 
net effect of the magnetic field is to reduce the local temperature distribu- 
tion. In the magneto-inviscid thermal layer the balance of thermal energy 
is controlled by the Joule dissipation and as this decays slowly to zero 
in consequence 7' tends slowly to T,. 

The magnetic field vectors are not without interest. In the magneto- 
viscous layer the axial component //, is for all practical purposes constant. 





ON AXIALLY SYMMETRIC STAGNATION POINT FLOW 401 


However, slow variations in both H, and H, occur in the magneto-inviscid 
layer. Using Tables 3 and 4 we find that to within a few per cent the 
components //, and H, are insensitive to m for 0 < m < 0-2. The dimen- 
sionless magnetic profiles h({), H(), and h’(f), H’(»), defined by equations 
(2.4) and (4.2), are shown in Fig. 4 by the representative curves for 
m = 0-1. 
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TABLE 4 
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SUMMARY 

The theory of steady axisymmetric solutions of the idealized hydromagnetic 
equations for a compressible gas is applied to an incompressible fluid. The aim of 
the paper is to investigate the solutions of the vorticity equation, which is a non- 
linear partial differential equation of the second order in y—the Stokes stream- 
function——and contains four arbitrary functions of ys. The solutions of this equation 
in some particular choice of arbitrary functions are obtained, among which is the 
generalization of Hill’s spherical vortex in the presence of a magnetic field. 

Generally both the velocity and the magnetic field vectors become infinite on 
a certain surface defined by the singularity of the vorticity equation. Section 7 
deals with this singularity. 


1. Introduction 

THis paper is mainly inspired by T. V. Davies’s paper (1) which was 
presented at ‘The British Theoretical and Applied Mechanics Colloquium’ 
held in Manchester, 6-9 April 1959. In that paper Davies deals with 
the case of a steady, compressible, inviscid fluid of infinite electrical 
conductivity, when both the magnetic and the velocity fields have an 
axial symmetry. Here that theory is applied to the case of an incom- 
pressible fluid. 

As the reduction of magnetohydrodynamical equations is nearly the 
same as in (1), except that here the density p is constant, the details of 
this reduction have not been repeated here. The aim of this paper is to 
obtain the solutions of the vorticity equation (2.10) in different cases and 
to describe each motion corresponding to these solutions. 


Generally both the velocity and the magnetic field vectors become in- 


finite on a certain surface defined by the singularity of the vorticity 
equation. Section 7 deals with this singularity. 


Chandrasekhar (3) has obtained the general equations, in terms of four 
arbitrary functions, as a particular case of his more general discussion. 


[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 4, 1960) 
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But there no reference is made to the singularity, nor to any general 
relation showing the dependence of all arbitrary functions on each other, 
such as the vorticity equation (2.10). 


2. The equations of motion 

The equations of motion of magnetohydrodynamics are the union of 
Maxwell’s electromagnetic equations and hydrodynamical equations, with 
the terms expressing the effect of interaction of fluid motion and magnetic 
field. In the case of a steady, incompressible, non-viscous fluid motion 
under the effect of a steady magnetic field these equations are 


] 
Vx(0xV) 4) jxH = v(2 + 4q?-4 0), vV.V=0, 
p p 
VxH = 4nj, V.j = 9, 
VxE=0, V.H=0, (2.3) 
j = o(E+V~H), (2.4) 
where V is the velocity vector, H the magnetic field vector, j the current 
density vector, & the gravitational potential, p the pressure, p the constant 
density, E the electric intensity vector, o the electrical conductivity, and 
q Vv. 
When the electrical conductivity is infinite (2.4) gives 
E —~VxH; 
hence (2.3) becomes Vx(Vx«H) = 0. 
Eliminating j between (2.1) and (2.2) we have 
] 
Vx(VxV) H (Vx H) = 0? 4 492+Q). 
4p 


Equations (2.6) and (2.7) are the basic equations to be solved under the 
conditions that the vectors V, H, and j are all solenoidal. In cylindrical 
polar coordinates r, 6, z, when V and H are azimuth-independent, it is 
possible to show that one can obtain a Bernoulli’s equation 


P y+ rH, f" 
p dnp 


and the vorticity equation 


F(), (2.8) 


ps Sa 
4npF’ 4 (049) ee) + a ED 
ve r “—< r2 < 


es ok’ +-vH,g'—rlf”, 
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or 


” j | ; ° ” 
brrpr? PY Al pp — hy tp + gg’ (yp? + 2) 
, 


, 


, K'(K — rgf’) +f" (gK —Anpr?f’) 4 AG rgf'\(gh —4npr?f \p 


2.10) 
where A g*?—4zp, & is the Stokes stream-function introduced by 
V.V — 0, fg, K, F are arbitrary functions of %, F’ the first derivative 
of F with respect toy, and &% isan unknown function of 7 and z. Alsou, v, w 
are the velocity components, //,, //,, 1/1, the components of the mag 
netic field Hin the principal directions in cylindrical coordinates, and 
suffixes denote partial derivatives of i. See (1). 

Now the vorticity equation (2.10) is a partial differential equation of 
the second order in ys, containing four arbitrary functions of 4, which are 
at our disposal. By choosing these functions suitably we can solve the 
equation and express ys in terms of r and z. Then 

=: r*gf Ws, 
r j*—4np r 
lyk — 4npr?f’ 


il, | 
. r g*—4np 


gy, 
* 





(rfl) , : (/1, /1,,), Js ae (r#1,) 


birp « 1 tarp cr 
2.11) 
follows. Bernoulli's equation (2.8) determines the pressure. 

It is clear that H- yV, if and only if 7, — gv, Le. iff’ — 0. Then the 
vectors V and H are everywhere parallel. The coefficient of the highest 
derivative in the vorticity equation is A = g?—4ap. On the singular 
surface g? 4p 0, which is generally a surface of revolution about the 
axis of z, both the velocity and magnetic fields tend to infinity. See 


section 7. 


3. Particular cases 

In this and the following sections we shall try to find particular solutions 
of the vorticity equation (2.10), which is a non-linear partial differential 
equation of the second order in y:. The equation contains four arbitrary 
functions of %, namely, f(y), gb), ACh), FQ), and &% is an unknown fune 
tion of r and 2 

Beginning with the simplest cases, we shall show that, choosing all or 
some of the arbitrary functions suitably the vorticity equation can be 
solved to determine ys. After determining the velocity field we shall give 
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a general picture of the motion and find the magnetic and the current 
density fields associated with the solution. 


3.1. Let K'=g’ =f" = F’ =0 


Let K = Ky, g = gy f' = fo. and F = BK, where Ky, go. fo, Ky are 
arbitrary constants. The equation (2.11) shows that v and //, are functions 
of r only. Now the vorticity equation (2.10) reduces to 


! 
Pre bet Pes = 9 (3.1) 


and we shall assume that g?—— 47p is not identically zero. 


(i) One well-known solution of this equation is (2) 


:\ 
yp = a (3 a) |: (3.2) 


t 


where a, ( are arbitrary constants; from (3.2) we can obtain the com- 
ponents wand w of the velocity outside the sphere r?} z* <= a*, which are 
finite at infinity. The component v of the velocity is, by (2.11), 

1 Ky r*gofo 

r ge Anp 
If we assume that fi, = 0, v becomes zero at infinity. Hence we have 


3 Uatrz 


Ky 
gi, -Anp r 
as 3 air? 


4 + 
(r?-4-22)' © 2 (r?4-22)!8 





wii 


The velocity is infinite at the origin, but if the origin is enclosed in any 
sphere of radius a, the velocity outside of the sphere is finite everywhere 
except on the axis. At infinity «= w— 0, and the motion which is 
irrotational everywhere reduces to a parallel flow in the direction of z, 
the magnitude of the velocity being U. 

Since / — 0 on the surface of the sphere, this surface is covered by 
stream lines. They are helices traced on the sphere. The stream lines outside 
the sphere are again helices, but as r increases both v and « diminish and 
the motion sufficiently far from the z-axis becomes parallel to the axis. 

Since f, ~ 9, (2.11) shows that H = g,V, that is the magnetic field 
vector is everywhere parallel to the velocity vector, and is a constant 
multiple of it in magnitude. The magnetic lines of force are the same as 
the stream lines whilst the current vector j is identically zero, 
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(ii) Suppose that y is of the form 
& = Rir).Z(z), 

R’—R'/r |Z" 

then (3.1) becomes — 1. 0. 
R Z 
The first term is a function of r and the second term is a function of z only; 
hence each term is a constant, m say, so that 
Z"—mZ : 
' _ (3.4) 
and R’— R'/r+-mR = 0) 
where m is an arbitrary constant which may be positive, negative, or zero. 
(a) If m n*, we have 

Z Ae™ + Be nz 
and the equation in & transforms to a Bessel’s equation by writing R = rP. 
where P is a function of r, namely, 

r?P"+rP’ + (n*r?—1)P = 0. 
Hence P Rir = Jnr), 
and R rJ\(nr). 
The second solution of Bessel’s equation is omitted since the velocity is 
finite on the axis. Hence 
wb rJ,(nr)( Ae" + Be-™) 


nJ,(nr)( Ae” ") | 
Yolo if K, 0 


Jo — 4p % 
uw nJ,(nr)( Ae nz} nz) 


n being chosen positive. The velocity is along the axis and is finite on 


and we have 


the axis, except at infinity; the arbitrary constants A, B, n, and m are 
determined by the boundary conditions. The motion is generally helical. 
(b) If m n®, we find 
Z = Acosnz+ Bsinnz, 
rI,(nr), 
and / rI,(nr)(A cos nz + Bsin nz). 
Hence 


nI,(nr)(A sin nz — Bos nz) | 
0 


Joto r mr if k, 


gi, —4np 
w = nl(nr)(A cos nz+- Bsin nz) 


where J, and /, are Bessel’s functions of imaginary argument. 
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(c) If m = 0, we find 


af = (Ar? B)(Cz+-D) ) 


C' 
so that u-— —--(Ar*+ B) 
r 


v=or ifk,=0 
w= 2A(Cz+D) 





where A, B, C, D are arbitrary constants and 
w = Jofo(gz—47p) = constant, 
which is also arbitrary. If B = 0 we find 


u= ACr 
v= wr , (3.8) 
2A(Cz+ D) 


When r = 0,u = v = 0, and w is only a function of z; hence the velocity 
on the axis is along it and a function of z. When z -DIC, w= 9, 
and the stream lines on this plane are equiangular spirals, w is positive on 
one side of the plane and negative on the other. The stream lines are 
given by the equation 

dr rdé dz 
ACr mr 2A(Cz4-D) 


The first two ratios give r = ae“ACi@e (3.9) 
the first and third give r(Cz+D) = b, (3.10) 


where a and 6 are arbitrary constants of integration; (3.9) represents a 
family of cylinders with their generators parallel to the axis of z, and 
whose projections on the planes z = constant are equiangular spirals 
given by (3.9) itself. The constant slope of these spirals to the radial 
vectors through the axis is —mw/AC. (3.10) represents a family of surfaces 
of revolution about the axis. The actual stream lines are the intersections 
of these two families of surfaces, the velocity being infinite at infinity. 
The magnetic field and electric current components are 


H,=—ACgr, H,= "Per, Hy = 2Ag,(Cz+D), 
Yo 
i 0, Je = O, js = Po, 
Jo 
where a = gofo/(g3—47p). 
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3.2. Let K'’ =y' =f" =0 and 4npF’ = (4np—gz)C, where C is an 


arbitrary constant 


The equation (2.10) becomes 
(3.11) 


One known solution of this linear equation is (2) 
dh ph, dal 'y?2(a2 -— r?2— 22), 
where a is an arbitrary constant. Since y — 0 on the surface of the sphere 
r? 2? — a*, the surface of the sphere is covered with stream lines: no flow 
passes the surface. But « and w become infinitely large at infinity, so that 
the solution is valid inside the sphere only and we take the steady flow 
viven by (3.3) outside. 
Now outside the region ys is given by 

(3.2) 
where 7? + 2? > a? and UU is the constant velocity parallel to Oz at infinity. 
This is a solution of the equation (3.11) when C = 0. The two expressions 
for & show that the normal velocity is zero on both sides of the surface 
of the sphere, and hence is continuous. The continuity of the tangential 


velocities require that C - 150'/2a*. Hence, inside the sphere, 


(3.12) 


where r* 2 


1 Ky r*gofo 
r qe tip 


In this case 


If Ay = fy = 0, then v = 0, and (3.12) gives the Hill's spherical vortex 
(2) inside the sphere r?4-2? = a®. And (3.2) represents an irrotational 
steady flow outside the sphere, which is finite and equal to U’ parallel to 
Oz at infinity. The motion is meridional everywhere. The stream surfaces 
are given by y — constant. Inside the sphere stream lines form closed 
curves in meridional planes around the points z = 0, r = a/y2 = 0-707a. 

Since fy = 0, we have H = g,V: thus the magnetic lines of force are 
coincident with stream lines. Current density lines are circles around the 
axis in planes z — constant, whose centres are on the axis. Its magnitude 
is — 150 gr S7a* inside, and zero outside the sphere. 

If Ay = 0, fy 4 0, then v = {go fy/(g2—4ap)}r = wr, say. This adds a 
rotation, with constant angular velocity aw, to the motion described above. 
The stream lines become helices traced on coaxial surfaces of revolution 
yy constant Inside and iby constant outside. 
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Magnetic lines of force are helices similar to the stream lines, but they 
are not coincident. 

If kK, 4 0. fo = 0, then v (1/r)K,/(g3 —42p). This adds a rotation 
to the meridional motion defined by ys, and y,, which is infinite on the 
axis and tends to zero when r + «. The stream lines are helices as before, 
the magnetic lines are the same as the stream lines since H = g)V. 

If one uses $4 

Yolo 
yy 47p 
Us . Ko outside the sphere, 
. r gi 4p 


r inside the sphere, 


so that v is finite everywhere, the v-components of the velocity and 
pressure on both sides of the surface of the sphere cannot be made equal 
everywhere, hence an elastic spherical membrane cannot remain spherical 
during the motion. 


4. Solution of the vorticity equation when 0 

Equations (2.11) show that in this case w — H, — 0. Hence the motion 
takes place in planes perpendicular to the axis of z and the magnetic field 
lines are in the same planes as the stream lines. 

When yf, 0 the equation (2.10) can be written 


, 


pee . = kK? 4 "e 
Sapr* kh } l(g? trp)? | ' ( 1) { tapr'(" 


2y( Kf’) 2K y'f'(g? + dnp) 
A A* ; 


Since F, f, g, and K are functions of y%, which is a function of z only, 


y2 


(4.1) 


(4.1) is satistied only if the coefficients of powers of r are the same on 
both sides. Thus we find 
yy" Big? — Amp), (4.2) 
. k? , 
(g? — 4p)yp? 4 C, (4.3) 


9 


g° —4mp 
tpg? —4np) F’ = g( Kf’) —Kf'g? 
Y 


24 4p 
q” dap’ 


(4.4) 
where B and C are arbitrary constants of integration and C 4 0. These 
three equations must be satisfied simultaneously. If F’ — 0, (4.4) is 
satistied by 
(a) f' = 4@, (b) g= 9, (c) K ; (4.5) 
and lastly if f’,g, K + 0 by 
(d) fgk = E(g (4.6) 


where £ is an arbitrary constant. 
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In case (a), F’=0, B=0 (ie. f’ = 0). 
The equation (4.2) gives f’ 0 and (4.4) is satisfied. (4.3) gives 


=[C(g? 4p) — K?} 
gq? ip : 


ws. (4.7) 


By choosing g(ys) and A(us) we can determine ¢(z), and hence wu, in terms 
of z. But since 


r g*—4A7p 


1 K | 
: ue 0 
r g?—Anp 


the stream lines are given by dr u rdé@v —dzw. Hence 


1 [ C(g?—4xp)— K?}} | 


(4.8) 


2 Ci, r Ge", (4.9) 
where k =| C(g?—42p)— K?|! K = a function of z only. Thus stream 
lines are spirals in the planes z = constant, which make an angle « with 
the radius vector at the point considered such that 

‘ rdé DK | 

an x a = = 
dr | C'(g°—4mp) -k?}! k 
is a function of z only. Thus this angle is constant in each plane perpen- 
dicular to the axis, but changes with z. The motion is real only if 

C(q?—4np)— K? > 9, 
which defines a region or regiou.s of validity. Since f’ 

im = gv. (4.10) 
i.e. the magnetic and the velocity field vectors are parallel; hence the 
magnetic field lines are coincident with the stream lines. The velocity is 
infinite on the axis and on the planes g?—4zp = 0. 


The eases (4), (c), and (d) give similar motions. 


5. Solution of the vorticity equation when yj. — 0 

The equations (2.11) show that u = H, = 0. Hence the stream lines, 
and the magnetic field lines lie on circular cylinders coaxial with the axis 
of z. They are generally helices, whose inclination to the generators of 
the evlinders is a function of r, and so constant on each cylinder. 

In this case the vorticity equation (2.10) can be written as 


Smpr? F’ +-r| (g? — 4arp)w? |’ 4 (=) + 4nprt(7) 


(5.1) 


[eae 2K f'g'(g?+ 47) 
A A? 
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where A = g*—4zp, supposed not identically zero, and w = (1/r)), is the 
velocity component in the direction of z. Dashes refer to differentiation 
with respect to ys, and & is a function of r only. Some particular cases are 
discussed in the next section. 

5.1. K = 0@ 

In this case (2.11) shows that 

v = (g 42p) My, 

i.e.v and H, are proportional, the coefticient of proportionality being a func- 


tion of r. The equation (5.1) becomes 


Sap F’ +-|(g?—4p)u?]'+ tpr( | = 0. (5.2) 


Now let us change all the derivatives with respect to ys, to derivatives 

with respect to r by using the operator 
d d 
is, : 

dus dr 

Hence multiplying by %, throughout and using the operator, we get 
. eal / fr 
Sap F.+-| w?(g?—42 -+ 4rrpr?{ —_-* __ - O. 
pF + [wig — tap), (a —4np)}, 

Now if we write u?(g?—4ap) = y, 
(5.3) becomes 


» 
yy? —Anpf;) + spy (f2)e— Set 2uF,| = 0. 
The formulae for v and H, are 


; rgf’ _ “en / | of, 
g’—4zp (g?—Anp)t, = (g? —4np)w 


H, — 4rprf’ 4rpf, 
GP Anp  (g?—4np)w 
Now equation (5.5) can be satisfied by 


» 


(> Si }- 2y F, 


| Y, =O 
and 
\ or y®—A4npf? 
(5.7) gives F, and from (5.8) we obtain 
La 
Oo = = —— 
| a(g? —4zp) |! 
F of, 


9= 


Ee 
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where g and f, are arbitrary, and a is an arbitrary constant. The stream 
lines are 
(5.11) 
that is, helices traced on coaxial cylinders r = c,, whose inclination to the 
generators v/w = gf, a is a function of r. 
The magnetic field and electric current components are 
An Eqa 
H, ’ H, ae ete \? A, = ——.—— } 
|a(g?—4zp) |! |a(g?— 4zp) |! P 
; z - (5.12) 
: ; HH is = —_(rH,) 
) ; 2 3r° a. rey 
1 2 io }3 jor - 


If we use (5.9) instead of (5.8) we get the set of equations 
0, y*—4npf? 0. 


The first defines F, the second gives 


y u?(g*—Amp) E. (4p) f,; 


hence , = (4zp)}| : tae } | 
g?—4np, 


q +f, }! 
(477p)'\g*—477p 


(5.13) 


a solution identical with (5.10) when ,(47p)f, = a. The stream lines 
(4p) - 
vere! r0+-¢, (5.14) 
q 
are again helices on the cylinder r c,. If for a certain value of r (say rp), 
f, = 0, and if [(g?—47p)/g*],_, 0. then both v and w are zero on the 
evlinder r = ry. If, for instance, 


\ | imp)a*(r -1o)*, 


and if we take + signs, we have 
uw \(4zp)a(r—To) | rages 
° (5.15) 
v = ga(r—r,) J 
Now if g(r,) is finite, both v and w become zero when r = 19; i.e. the 
velocity is zero on the cylinder r = ry, and has opposite signs inside and 
outside this cylinder. 
The stream lines are given by (5.14) everywhere. If we assume g(r) 
and a < 0, the velocity on the axis is 


w —,(47p)aro, 
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i.e. the fluid particles on the axis have only a translational motion along 
the axis, which is constant and positive. As r increases, w decreases and 
tends to zero as r->ry. On the other hand, since g(r) = r, and a < 0, 
v increases first, takes a maximum value, and then decreases again to zero 
as rfp. 

Outside the cylinder, r > ry, both v and w change sign to become 
negative. The velocity is zero on the cylinder r = rg, but for r > ro, v and 
uw become negative and increase indefinitely as r—> x. 

The components of the magnetic field and electric current are 

H, = 0, H, = 4xpa(r—ry). H, = ,(4zp)ar(r—ry) 
; p\h ; 2r—ro . (6.16) 
ji = %, a= ar, Js = pa 

(2) r 

Hence the magnetic field vector H is also zero when r= ry. The 
current vector is infinite when r = 0, it is constant when r = ry, and j, 
changes sign when r = 4}7r,. 

Again, if we assume that 


, (4zp)a? sin?r, 


and if g() is finite, by taking the + signs only we have (see (5.13)) 

w = ,(4zp)asinr | ees 

; (5.17) 

v = gasinr | 

The velocity is zero when r = 0, 7, 27, 37,..... Hence the cylinders r = kz 

(k integral) can be taken as boundaries and the motion between them is 
then known. 

The stream lines are everywhere the helices given by (5.14), on cylinders 
coaxial with the axis of z. Their inclinations to the generators of the 
cylinders are constant and equal to vw = g/,(47p). The total velocity q¢ 
is given by q = (v?+-w*)! = (g?+ 4xp)tasinr. 

Therefore the velocity is zero when r = 0, 7, 27,... and a maximum 
when r= 7 2, 37/2, 57/2,..... The maximum value of the velocity is 
(g?+-47p)'a, which is a function of r provided that g is also. The velocity 
changes its sign as each bounding cylinder (the cylinder on which the 
velocity is zero) is passed. The magnetic field and current components are 


H, = 0, H, = 4zpasinr, H, = \(4zp)gasinr | 


, gece. p\h ; , p~ » (5.18) 
i). o=— — ga cosr, == " + Ss 
Ji J2 (; ga cos? Js pay ; cosr) | 


7 


The magnetic field vector H is zero on the bounding cylinders formed 
of the stream lines, and a maximum on the intermediate cylinders. 
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§.2.9=0 

In this case g?—4zp is always different from zero and there is no 
singularity. Then (2.11) shows that H, = H, = 0, H, = —rf’. Therefore 
magnetic field lines are circles about the axis of z, instead of helices. 

The equation (5.1) becomes 

Ald 9 avr l ra\r rayr 
Sapr?F’ —4rpr2(w?)’ = r [(K?) +4apr4( f’*) |. 
7p 


and after transforming all the derivatives from u% to r as in section 5.1 


this can be written as 
(4rpr)?(2 F —w?), = (K?),+- txpril fr :) ; (5.19) 
rue), 

; :& : ~ 9 
If constant, (9.20) 
ru 
then (5.19) becomes 

(4rpr)?(2F —w?), 


and (5.20) gives 


whence, from (2.11), 
dopr 


f.and K are arbitrary functions of r and (5.21) defines F. 
If for example we choose F and A so that F = bf and K 


relation (5.21) becomes 


ofan, —2hle 202] — ag, 


or 6 @eF 
Hence either f, 
= a*br?. (5.23) 
if f. = 0, we obtain w = 0, v = (1 r)f,, which represents a pure rotation. 
If f, ~ 0, we have to solve (5.23). If we put f = rP, this becomes 
r . 
rPP+rP.+ (ar?—1)P = abr’, 


which is Bessel’s equation of order one. Thus 


P , cJ,(ar)+br, 


r 


where a, b. ¢ are all constants and 
» 


2h 
w cJ,(ar) + — | 
a 


cJ, (ar) ihr 
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|The second solution of Bessel’s equation is rejected since w is finite on 
the axis.| The magnetic field and current components are 


H, = 0, H, = - 
(5.25) 


ji = 9, je = 9, a= on 


Both the magnetic field and current field are very simple. Magnetic 
field lines are circles with their centres on the axis in planes perpendicular 
to the axis. The current vector is constant and parallel to the axis. The 
solution (5.24) can be obtained directly from (5.21) by using F’ = 6, 
f" = 0, K = 4zpay. 

$.2. f' = 0 

The equation (2.10) shows that H, = gv: and since H, = gu, H, = gw, 
the magnetic field vector is parallel to the velocity vector everywhere. 
Therefore the stream lines and magnetic field lines are coincident. 

Equation (5.1) becomes 


iia 
Sapr?F’ + r?{ w?(qg?— 4p) |' 4 == @, 5.26 
p [u2(g2—4p)] L : ‘as) (5 


d d we get 
‘di dr © 


or if we transform the derivatives from y to r by & 


SapF,+-|w?(g?— 4p) |, + ; = = 0, (5.27 
' - \g? —4p r 


where w = (1 r)s, is the velocity component in the direction of z. 


Also v : --{- x \ 
r\g-—Azp, 


and ._=- = H, = gu, 


l nit oe gk 
=” J3 ~ Gor g*—47p/], 





Hence, as was stated at the end of section 2, 
H = gV, (5.29) 
where g is a function of r. 
The equation (5.27) is a relation between four arbitrary functions F, g, 
K, and w. It gives w in terms of the other three functions in the most 
general way, namely, 


u?(g?— 4p) = C— Sap F — | (a5) a: (5.30) 
. r 


where C is a constant of integration, F, A, and g are arbitrary functions 
of r. 
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kK2 
Let - ar™, 
g°—4rp 
where a and m are arbitrary constants. Then (5.30) gives 
2 Sap F —mar™-?/(m—2)\} > 
g°—4mp | 


. ( a }! 
and ) fim —2)/2) _— : 
g° — 47p; 


In this solution F and g are arbitrary functions of r and a, m, C are 


(5.31) 


arbitrary constants. The stream lines are again helices traced on coaxial 
circular cylinders r = ¢, about the axis. Their inclination to the generators 
and their pitch are generally functions of r. 

If, for instance, m 4, F = 0, (5.31) becomes 


'— 2ar?\} - a }! 
= “* v= Frl\— “ 
q- ae 9° — 4zzp, 


Ca must be positive. Generally this defines a motion inside and on the 
surface of a certain cylinder r < (C/2a)'. Also 


H, = 0, H, = gv, H, = gw 


(qu) : (rgv) ; 
; a — qu)... e=- qv), 
Ji Je * gw), J3 4 q | 


4 or 


a 


6. Motions in meridional planes 
These are obtained by assuming v = 0. But since 
1 K—?rof 

i; gq? 4p’ 


we have A—r’gf’ — 0. Hence either 
(i) K q o. or 
(i) A ow f* = 8. 


In the first case the equation (2.10) becomes 


l nai » 
$7 Ur hb ta.) Lif’ f 4rpr*F’, 
: 2 2 


and in the second it becomes 


> : | , 2; 9 on ai 
(g* 47p)( b+ b..] L-gq'(p2+p?) = - 4zpr°F ; (6.2) 
\ r = ® 
Solutions of (6.1) 


When f” = F’ = 0, (6.1) reduces to (3.1). Hence all the solutions of 
(3.1) obtained in section 3, with the additional condition A = g = 0, are 
also solutions of (6.1). For example (3.5), (3.6), and (3.8) when g, = 0. 
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When f” = 0, F’ = C, (6,1) reduces to (3.11). Hence (3.12) is a solution 


of (6.1) when g, = 0. 
Also all the solutions obtained in sections 4 and 5 represent meridional 


motions when v = 0. 
Solutions of (6.2) 
Let & = R(r).Z(z), g?@—42p = 2, ie. gg’ = 1, and F’ = 0, Then (6.2) 


becomes l 
2Rz( RZ — RBL+ rz’) (R'Z)2+-(RZ'P = 0. 


Dividing by R?Z? we find 
2Z" Z') 2k" 2 R%_ 
| rt Ror ztB) gies 


The first group of terms is a function of z, and the second is a function 


of r only. Hence each must be a constant, m say, whence 
(6.3) 


2Z2Z"+Z"—mZ? = 0 
9 

2RR’—=— RR'+ R? +m = 0 
r 

are easily found to be 


rhe solutions of these equations, when m = 0, 
Z = (Cz+D)i, 
R = (Ar?+ B)i; 
yb = (Ar?+ B)i(Cz+ Di, 


2C(Ar?+ B)i 


4 = ~ 3r(Oz4+-D)i 
_ 4A(Cz+ D)i | 


= 3(4r24+- By 


hence 
where A, B, C, and D are constants. Thus 


0 = 


The stream lines are given by y% = constant, i.e. 
(Ar?+ B)(Cz+-D) = constant. 


Compare this with section 3.1 (ii) (c), i.e. with the expression for % in 


(3.7). The stream lines are the same in both cases. 
The velocity is infinite on the axis, and cn the surfaces 4 = 0, owing 


to the singularity of the equation, 
H, = (2y+-4np)!u, H, = 9, H, = (2b+-47p)'w 
: F ] 
A= 0, a= 4 (H,,—H,,), 
7 


j, = 0 
showing that the magnetic field is poloidal, but the current density field 


is toroidal. 
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Since the solution of the equations (6.3) when m = 0 is the {rds power of 
the solution of the equation (3.4) in the same case (i.e. when m = 0), one 
may expect that this relationship between the two sets of equations exists 
when m + 0. In fact this is true if we assume m = $n? when m > 0, and 
m —{n* when m < 0 in (6.3). 
Hence the solution of (6.3) when m = §n? is 
Z (Ae nz1 Be- nz)§ 
R = [rJ,(nr) }*; 
whence ys = [rJ,(nr)(Ae™+ Be-™) }f, 
and the components of velocity are 


__ 2n{J,(nr)|*( Ae — Be-"*) 


3rh(Aen* | Be nz) 


2nJ,(nr)( Aen + Be-)i 
: 3r*| J,(mr) | 


(6.6) 


where J,, J, are Bessel’s functions of order one and zero respectively. 
The solution when m jn? is 


Z = (Acosnz+ Bsinnz)i, 
R |rZ,(nr) }§, 
so that ys = [r],(nr)(A cos nz+- Bsin nz)}F. 


Hence 2n| 1,(nr)|*(A sin nz— B cos nz) | 


3r4(4 cos nz+ Bsin nz) 


(6.7) 
2nI,(nr)(A cos nz+ Bsin nz)! 
3r4{ L,(nr) |! 


where /,, /, are Bessel’s functions of imaginary argument. In all cases 
(6.4), (6.6). and (6.7) the stream lines are given by y% = constant. They 
are the same as in section 3.1 (ii). The magnetic and current density 
fields are found by (6.5). 


7. A discussion of the singularity of the vorticity equation (2.10) 
The coefficient of the highest derivative in the partial differential 
equation (2.10) is. A = g?—4ap. Hence, when g?—47p = 0 we may expect 
an infinity in the solution. In fact (2.11) shows that as g?—47p > 0 both 
v and H, tend to infinity. 
Consider the equations giving expressions for v and H,, namely, 
4rprv—grH, = kK, (7.1) 
grv—rH, = r°f’. (7.2) 
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If g?—4zp ~ 0, these define rv and rH, as 


rH, = gk —— 
g — 7p 


If g?—47zp = 0 we cannot sotve for rv and rH,, and the two equations 
(7.1) and (7.2) are only consistent if 
K—r*gf’ = 0; (7.5) 
the two equations then become identical. Instead of the two relations, 
we can take the single relation (7.2), which can be written as 
gvu—H, = rf’. (7.6) 
The Bernoulli's equation (2.8) does not change, and the vorticity 


equation, if we take the form (2.9), becomes 


4apF’ +9 (y2+-2) = ok 4 vH,g' —rH,f’. (7.7) 


2 
Using g?—4zp = 0 and (7.5) it follows that 
gk —4Anpr2f’ = 0; 

hence both v and H, become indeterminate. We shall separate the two 
cases g’ = Oandg’ + 0. 

7.1. 97 =0 

This means that g is an absolute constant, and the relation 

g?—Anp = 0 (7.8) 


is satisfied whatever the values of r and z, i.e. everywhere. In this case 
the vorticity equation (7.7) becomes 


4npF’ = !vk’—rH,f’. (7.9) 
r 


The equations (7.6) and (7.9) are two equations to determine v and H,. 
If K’—r’gf” #0, v and H, are 
4apF’ oT. | 
KR igf 
. ; (7.10) 
H 4rpg F’ fh’ ‘ 
= kK’ rgf” 
Equation (7.5) shows that A and f’ are either both zero or both different 
from zero. If K = f’ = 0, (7.10) does not define v and H,, since 


K’—r*gf” = 0. 
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But since f’ = 0, (7.6) shows that gv—H, = 0, that is H, = ,/(47p)v by 
(7.8). Hence it follows from (2.11) that 
H — , (4zp)V. (7.11) 
This is the relation given by Chandrasekhar (3) as the simplest solution 
of the magnetohydrodynamic equations under the conditions considered. 
The equation (7.9) then gives F = constant; hence Bernoulli’s equation 
° . 
(2.38) 1S p 
Pp 


+-4q?+Q = constant. (7.12) 


The results (7.11) and (7.12) follow directly from the main equations (2.6) 
and (2.7). 

On the other hand, ifone or other of A’ or f” is not zero, (7.5) indicates that 
ws is a function of r? only. Then u = H, = 0, w = (1/r)b,, Hy = (1/r)go¥,, 
v, and H, are found by (7.10), where gy = ,(47p). 

If forexample A = gyW?. f’ = ab, 4apF’ = aby, where a, b are arbitrary 
constants: from (7.5) we get 

Jo vs(us —ar?) 0. 
Hence either 0. or &b ar, 

We cannot accept the case & = 0, since then K = f’ = 0 and (7.10) does 

not define v and H,. When y% = ar? we have 
u 0. ’ elf ar), == Za, (7.13) 


Yo 
The stream lines are then 


2, (4727p) , 
1» . 2 iy T 
»—ar- 


ce Co, 


i.e. helices on coaxial cylinders, whose inclination to the generators of the 
evlinder is a function of r. Also 


H, 0), H, = ar(b—2ar?), 


A=9% =O j= — 


Bernoulli's equation gives p as a function of r. 
1a o£ +8 


If g is a function of ys, g’ is not generally zero, and the equation 
g?—4np = 0 is satisfied only on a surface of revolution defined by itself. 
Outside cf this surface g?—4zp + 0, and the normal equations (2.10) and 
(2.11) are used. On the singular surface g?—4zp = 0, we have to use (7.5), 
(7.6), and (7.7) instead. 

If A and f’ are identically zero, v and H, will be zero everywhere. We 
can assume that this is true even on the surface of g?—42p = 0. 
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If K—r’gf’ is not identically zero, it must contain g?—47p as a factor, 
since it must vanish when g?—47p = 0. Hence 
K—r*gf’ = a multiple of (g?—4zp)", 


and thence K = a(g2—47p)” ) 


f’ = Bg?—4np)™ J 
where « and § are arbitrary functions of y% which are finite when 
g?—An 0, 
and m is a positive arbitrary constant to be determined suitably. 

This fact is realized in the general equations of sections 4 and 5. For 
example, the equation (4.1) suggests that both A? and f’* must be multiples 
of g?— 4p. It is seen in (4.2) that iff’ + 0, it must be a constant multiple 
of g?—4rp. (4.3) shows that if A = 0, it is a multiple of g?—4zp = 0. 
Again, the equations (5.1), (5.2), and (5.27) show that A and f’ are either 
zero or multiples of g?— 4p. 

In any problem we usually know g?—4zp as a function of ys, also we 
know the value of & in terms of r and z, i.e. the singular surface g?— 4p = 0 
is known. If on the singular surface the formulae (7.3) and (7.4) become 
indeterminate we use (7.6) and (7.7) to determine v and H,, under the 
conditions (7.15) imposed on K and f’. 

Usually, if A and f’ have the form (7.15) the formulae (7.3) and (7.4) will 
be determinate and can still be used to determine the values of v and //,. 


But in general the velocity is infinite on the singular surface. For example 
take the solution (4.8). If A = a(g?—4zp), then f’ = 0, and v will be 
finite on the singular surface but w is infinite. In (5.10), if f, is chosen 
to be a multiple of (g?—47p)!, K = 0, v becomes finite, but w = 0 when 
g?>—4zp = 0. On the other hand, in (5.13) it is possible to choose f, as 
a multiple of g?—47p, so that both wand w become finite on the singular 
surface, and everywhere if g is finite everywhere. 
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SUMMARY 

A theoretical analysis is made of the flow of an incompressible non-Newtonian 
viscous liquid in a curved pipe, looking for differences in observable characteristics 
from the corresponding case of Newtonian flow. Such a motion is of interest to the 
experimentalist because the flow could be readily attained and controlled in practice, 
the most easily measurable quantities being the axial pressure gradient and the 
volume rate of flow. It is assumed, for the purpose of mathematical analysis, that 
the curvature of the pipe is small, more precisely that the radius of the circle in 
which the central line of the pipe is coiled is large in comparison with the radius 
of the cross-section, A solution is developed by successive approximations, the first 
approximation corresponding to the flow of a Newtonian viscous liquid as given by 
Dean (1). The streamlines in the plane of symmetry and the projection of the 


streamlines on a normal section are compared with those of a Newtonian liquid. 


1. Introduction 

LIVLIN (2), (3) has proposed that the macroscopic behaviour of certain 
fairly dilute solutions of high polymers can be represented by equations 
of state in the formt 


Ss »  g 
Pik PO ip + 26 jp + Leis CRs (1) 
where the viscosity coefficient » and the normal-stress coefficient ¢ are 
scalar functions of the flow invariants 
i, », 1, = §(e,)?—enen|, Ig = det(e,,). (2) 


Here p,. and e,, are the stress and rate-of-strain tensors, P is an un- 
determined isotropic pressure which can be superposed on an element of 
liquid (in any state of motion) without affecting the flow, and 4,, is the 
Kronecker delta (8,, lift —kandd, = 0ifi #~k). A formula of the 
form (1) in which the coefficients » and ¢ are physical constants has been 


discussed by Reiner (4); Reiner considers the material to be compressible 


+ The usual convention of summation applies to repeated small Latin suffixes except r 
and z. 


{Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 4, 1960] 
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and accordingly the quantity P is a function of the density, i.e. of the 
dilatation. 

A variety of interesting theoretical predictions concerning liquids charac- 
terized by the rheological equations of state (1) are revealed in some 
simple types of laminar flow discussed by Rivlin (2), (5) and by Braun 
and Reiner (6), results which are not a feature of the quasi-linear theory 
(¢ = 0). In steady rectilinear flow in a straight pipe of circular cross- 
section the incompressibility condition J, = 0 is consistent with a velocity 
distribution [0,0,W(r)] when referred to cylindrical polar coordinates 
r, d, in which the z-axis coincides with the axis of the pipe. The equations 
of motion require (5) 
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where » and ¢ are in general functions of r, and A and B are arbitrary 
constants. The non-vanishing physical components of stress are then 
found to be 
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It is seen from (3) and (4) that the axial gradient of 22 is a constant A 
and that the differential equation for W(r) is independent of the normal- 
stress coefficient C. 

For reasons of mathematical convenience we shall assume throughout 
the ensuing discussion that the viscosity coefficient » and the normal- 
stress coefficient ¢ are physical constants. Equations (3) reduce, in this 
case, to 1 342 
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where we have taken r = a to be the boundary of the pipe. 


2. Flow in a curved pipe 

We consider the steady motion of the liquids characterized by (1) 
through a pipe of circular cross-section, of radius a, with the line of centres 
coiled in a circle of radius 6. The coordinate system used is that adopted 
by Dean (1) in discussing the flow of Newtonian liquids under the same 
conditions and is illustrated in Fig. 1. Ov is the axis of the circle in which 
the pipe is coiled, C is the centre of the section of the pipe by 2 plane 
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through Ov that makes an angle @ with a fixed axial plane, and CO is 
perpendicular to Ov and of length 6. Any point P of the section 6 = con- 
stant is referred to by the orthogonal curvilinear coordinates r, ¢, 6, where 








\ 
ur, 


Pic. 1. The coordinate system (r, d, @) chosen to describe motion in a curved 


pipe of circular cross-section. 


r is the distance CP and ¢ is the angle CP makes with the line through C 
parallel to Ov; the line element is then given by 


ds* dr? -r* dd? +-(b- rsind)*dé?. (6) 


It is assumed that the motion is steady, and, further, that the velocity 
components in the coordinate directions r, ¢, 6 may be written U(r,d), 
Vi(r,d), W(r,d), independent of 6. The rate-of-strain tensor has then the 
physical components 
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and, from (1), the stress tensor has the physical components 
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The stress equations of motion require, in the absence of body forces, 
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The equation of continuity, J, = 0, requires 


cu U Usind ] oV V cosd 
cr r 'b+rsind | red | b4 rsind 


(12) 


If the pipe were straight a/b would vanish, and equations (8)—(12) would 
be satisfied by the expressions (5) for W and P together with U = 0, 
V = 0. 

We now assume that the curvature of the pipe is small, that is, a/b is 
small. Following Dean (1) we write 
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where u, v, w, and p are all taken to be small, of the order of ab, and 
u,v, ware independent of 6. Neglecting terms of order (a/b)? the equation 
of continuity becomes 
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and is the condition for the existence of a stream function y% such that 
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The stress equations of motion (9)—(11) now reduce to 
— — (a? sin d— a y, (30° ~Tr?)sind 
167°) 
op 1 wt i 1 ep) CA/, ew y 
on Y + +- o = 
cr — ror rod? cr? 
pA? "A 
16," )? cos d6— 37% 7r?)cos d 
_leép Ley | l a) tA ¢ ("" 
Ly 
or 


r&b or* ror "72 ed?  -2n &d ér * 


rsing 


27 


pA cu 3A 
2n ch 2b 
cp ew lew 1 @w\ CA 6 (eb led, 1 ey 
ae eee ae tie tama =«((18) 
or? ror r* ad? 2y od\or? rer | r= ed® 
where we have written, in conformity with (13), z = 60 and (1 b)¢/c@ = @/ez. 
It will be seen that there exists a solution of equations (16), (17), and 
(18) in the form 
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where ws’, w’, and p’ are functions of r/a only. sniiaclieas for %, w, and 
p in terms of the dimensionless quantities ys’, w’, and p’ in (16)-(18), the 
following equations are obtained: 
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where # (a Reynolds number) and S are the dimensionless parameters 
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A is the operator defined by 
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and 7’ is written for r/a. The condition that there is no slip at the boundary 
of the pipe is satisfied if 
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Also w'/r’, di’ /dr’ and w’ must be finite at all points in the region 
0O<rcl. 

Dean (1) has shown, in the case when S = 0, that a necessary condition 

for a solution of physical significance is 


i = dy at r’ = 0 
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and w' =p =0 atr’ = 0, 
and his argument applies generally when S + 0. 
Eliminating p’ from equations (20) and (21) we obtain 
A(Ay’ +-2Sw’) = 4R(1—r’2)r’+28Sr’. 
Comparing equations (22) and (28) it is seen that w’ may be eliminated 
quite easily; we find that 
A*p’ -4S2Agb’—4RSp’ = 4Rr'(1—r’?)4+ 1687’, 

and this equation for y’ is to be solved subject to the boundary conditions 
25) and (26). When S is sufficiently small a solution of (25), (26), and 
(29) can be developed with successive approximations for ys’ equal to the 
partial sums of a series 

fs’ thy ' Sys, t S7ubs, ' SB, eo (30) 
in the following way. The series (30) is substituted in (29), and the left- 
hand side is expanded in ascending powers of S. Equating in turn the 
coefficients of S°, S', S?,... differential equations for Yo, ¥,, Yo,... are 
obtained. The process of substitution and equating coefficients yields the 
following equations: 
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Each successive approximation to us’ must satisfy the boundary condi- 
tions, and this requires 
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also y,,/r’ and dy, /dr’ (all n) must remain finite in the region 0 < r’ 
The first approximation yf, is found to be 
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in agreement with Dean (1). Substituting this value of %» in equation 
(32) we obtain R 
Ath, = r’'(1—r’?)?(4—r’?) 75+ 16r’, (37) 
the solution of which is 
R2 
4+ br’s 
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where A, L, M, and N are constants. The boundary conditions (34) and 
(35) require L—=M—O0 
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The resulting expression for yf, is 
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There is no difficulty, in principle, in finding the remaining coefficients 
i, wy, ete., but the algebra becomes rather heavy. For the purpose of 
illustration we shall now assume that to a sufficient approximation ’ 
may be taken as a linear function in S, 
By inspection of equations (20) to (22) it follows that we may write 

Pp’ = Pot Spy +S*pyt+ S8pg+... (41) 
and w’ = wet Sw, + S2w,+ S8ws+.... (42) 
The boundary conditions (25) and (27) require 
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n 
also w,, (all x) must remain finite in the region 0 < r’ < 1. Substituting 
(30) and (42) in the differential equation (22), and equating in turn the 
coefficients of S” (x 0,1, 2,...) on both sides of the equation, we obtain 


Aw, = 6r’—2 Ry, (45) 
Aw - —2Ry,, t 2Ay,,, 1 (n > 1), (46) 


where yo, Y,... are already determined. Substituting for %, and yf, in 
(45) and (46) the first two coefficients of the series (42) satisfy 
R2 
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and 
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The solution of (47) satisfying the conditions (43) and (44) has been found 
by Dean (1) to be 
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Integrating (48) under the conditions (43) and (44) we obtain 
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Substituting the series (30), (41), and (42) into (21) and equating in 
turn the coefficients S" (n — 0,1, 2,...) we find 
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The first two coefficients are found to be 


R 


J! 872 | Dyld 
Po r (9—6r'?+4-27 Mo” 


2 
yy by’ (55 —93r"2) 4-9’ (—212 4+ 6472 — 72r’4 4 32r’6— 4r’8) ce (55) 
There is no difficulty in principle in proceeding to evaluate the remaining 
coeflicients in the expansions (41) and (42), but it will suffice for illustra- 
tion purposes to neglect terms of order S?. 
What is perhaps of most interest is the effect of the normal-stress 
coefficient € on the relation between the axial pressure gradient and the 
rate of outflow. The rate of outflow is 


| a®r’W dr'dd 
r<0 df 
and, to the above order of approximations, is independent of w’ (and 
therefore of ¢), so that the relation between axial pressure gradient and 
rate of outflow is the same as if the pipe were straight. Dean (7), in a later 
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paper, has shown that this relation will be modified if terms O(a/6)? are 
included. Thus, to determine how the rate of flow depends on the normal 
stress coefficient ¢, for a given axial pressure gradient, we should have to 
consider the terms O(a/b)*; this has so far proved too difficult. 

The flow lines will vary with the dimensionless parameter S in a way 
that can readily be demonstrated even when terms O(a/b)? are neglected. 
The differential equation for the streamline in the central plane is 

or = Crane’ = +4), (56) 


and to a sufficient approximation we may write (56) as 
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the sign + being the sign of ¢. It follows, on integration of (57), that the 
equations of the streamlines are 
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6 being measured from the point where the streamline crosses the central 
line (r’ -- 0). The value of @, is in agreement with that of Dean (1), 
corresponding as it does to the value of @ when S — 0 (a Newtonian 
liquid). The (@,7r’) relation is seen to be independent of the ratio a/b. 

For a given value of 7’, 6 varies with the dimensionless parameters R 
and S; in the case of a Newtonian liquid (S 0) 6 varies inversely as R, 
The value of 4 increases steadily with r’ and tends to infinity as r’ tends 
to unity, while @, decreases steadily and tends to minus infinity as r’ tends 
to unity. However, in the range of values of S for which the approxima- 
tions made in this paper are valid, the function 6,-+-S@, tends to infinity 
as r’ tends to unity; the speed at which this function approaches infinity 
increases as S decreases. This shows, quite generally, that the curvature 
of the streamlines in the central plane increases (so that they follow more 
nearly the line of the pipe) as S decreases. Numerical illustrations are 
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now given for the particular boundary and Reynolds number considered 
by Dean (1), namely, 


R 63-3, ab ] (61) 


3? 


and for different values of the parameter S. (As Dean (1) has pointed out, 
the value one-third for the ratio a} is rather large and it is doubtful 
whether the approximations made in the above discussion are strictly 
valid for this value of a/b.) The values @, and @, for corresponding values 
of r’ are given in the Table. It is seen that the contribution of 6, to @ is 
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positive when S < 0 and negative when S > 0. There is a restriction on 
the magnitude of S, if the expressions for ’, w’, and p’ are to be con- 
vergent. [t is assumed that the series (30), (41), and (42) are all convergent 
when S is less than about 1-5. as seems likely from the first few terms. 
The forms of the whole streamlines for the particular values S —1-0, 
S = 0, and S 1 are shown in Figs. 2, 3, and 4, respectively. The 
effect of a positive S (a liquid which can exhibit the positive Weissenberg 
effect) is to decrease the angular distance that the fluid particles in the 
central plane travel in going from points near the inner edge of the pipe 
to points near the outer edge: the effect of a negative S (a liquid capable 
of exhibiting the less frequently observed negative Weissenberg effect) is 
to increase this angular distance. The values of S have been chosen rather 
large in order to demonstrate the qualitative effect of non-Newtonian 
behaviour: in an accurate solution we should, of course, have to include 
more terms in the series expansion. 

It is of interest to draw also the curves of intersection of the surfaces 
uw — constant with a normal section 6 = constant. Dean (1) has done this 
in the case of a Newtonian liquid (S = 0). The curves have the polar 
equation 


sec d kr’(1—r’?)?(1— }r’?) » 
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(1—}r’?) ‘i 48r L3r "5600 Ry (62) 


where & is an arbitrary constant. When S = 0 the solution is in agree- 
ment with Dean (1) and the relation between r’ and ¢ in this case depends 
neither ona} nor on R. For the class of liquids considered here it is seen 
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Fic. 2. Streamline in the central plane when S —10and R = 63:3. 


Fic. 3. Streamline in the central plane when S = 0 and R = 63-3. 
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that for any given value of r’ the value of sec ¢ varies with the Reynolds 
number FR and with the parameter S; again the (r,4) relation is inde- 
pendent of the ratio a/b. We shall take the Reynolds number to be fixed 
and to have the value 63-3, and consider the effect of a variation in S, 


Fic, 4. Streamline in the central plane when S LOand R — 63-3. 


From equations (13), (15), (19), (80), (36), and (40) we see that to a 


sufficient approximation V vanishes, for all values of 4d, when 
4—23r'2+-7r’s 


S| (91—599r’? + 401 r’4— 1247'S + 1 17’) R 1 (1 br’2) = 0. (63) 
2400 R 

In the case of a Newtonian liquid (S 0) the only relevant solution is 
r’ =~ 043 (in agreement with Dean (1)) independently of the Reynolds 
number F&. Taking R 63-3 the relevant solutions of (63), in the parti- 
cular cases S l-4 and S 1-4, are r’ 0-58 and 7’ 0-42, respec- 
tively. Thus at the points r’ — 0-58, 6 = nz when S l-4, 7° = 0-43, 
&b nz when S 0 and 7’ 0-42, 4 naz when S 1-4, both Ul and 
vanish, and there is a particular streamline in the form of a circle 

y constant, 4 nz. For any particular value of S, there is therefore 
a limiting surface ys, (the constant 4, dependent on S) which takes the 
degenerate form of a single circular streamline in a plane parallel to the 
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central plane. The intersection of the streamlines y = ys, with a section 
@ = constant are denoted by dots in Figs. 5, 6, and 7. The line % = y, is 
defined by k = 49-8 when S = —1-4, by k = 3-7 when S = 0 and by 
k = 1-9 when S — 1:4; k = m, for all values of S, corresponds to the pipe 
wall. Denoting the value of k corresponding to ys = y, by k,, the surfaces 
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Fic. 5, Curves of intersection of the surfaces ys = constant with a normal 
section # — constant when S l-4and R — 63-3. 


_— 





Fic. 6. Curves of intersection of the surfaces ys — constant with a normal 
section 6 — constant when S — 0. 


corresponding to kk, = 1, 1-1, 1-3, 1-8, 3-5 and o (along the pipe wall) 
are indicated in Figs. 5, 6, and 7. It is seen that as S increases from 
—1-4 to 1-4 the distance from the central plane of the limiting streamline 
ys — w, decreases from r’ = 0-58 to r’ = 0-42, an effect which might be 


clearly observable in an experiment; the general tendency is for the 
surfaces ys — constant to crowd nearer the limiting circular streamline as 


S decreases. 
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We have already shown that the angular distance (@), in which the 
projection of the fluid element describes part of the closed curves such as 
those in Figs. 5, 6, and 7, decreases with the parameter S for a fixed 


>. 
a 





Fic. 7. Curves of intersection of the surfaces Ws constant with a normal 


section 6 constant when S l-t4andR 63-3. Broken lines indicate 


the corresponding curves for a Newtonian liquid (S 0). 
Reynolds number R. This may be shown also from the (¢,@) relation. 
To a sufficient approximation we find that 
dé PSS?" 2PRSNr’ 
cos ¢ o = , io - ’ T ” m 
dd R(4—23r'2+7r'4) R(4—23r'2+7r’!)? 
‘is a 24 
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| 2400 R 
where the relation between r’ and ¢ is known. For points very near the 
boundary we can write r’ 1 approximately and (64) becomes 
dé 
dd 


(1 Sy'8) (64) 


0-38 + 0-23S)see d. (65) 


The angular distance in which a fluid particle near the boundary goes 
from a point where ¢ 1 to d = 0 is given by 
0 
| (—0-38-+-0-238)see ¢ dé = (0-38—0-23S)log tan(}z 
. 
This shows plainly that [a\? ” increases as the parameter S decreases, 
and tends to infinity, whatever the value of S in the range S 1-5, say, 
as a tends to $7; these results are compatible with those derived in the 
earlier part of the discussion. 

It is of interest to pursue the problem further to illustrate the kinds 
of stress distributions which arise in non-Newtonian flow. The stresses 
acting across the bounding surface (r’ = 1) are |—rr],_, and [72],_,. 
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Neglecting terms of order (a/b)? we obtain 
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where 77, and 72, correspond to the values of 77 and 72, respectively, when 
a » vanishes (see equations (4)) and are functions of r and z only. If 6D 
is the axial drag on a small axial length 5z of the cylinder then 
o=27 
8D = | Radddz = 2naFzy dz, (69) 
=o 
exactly as for a straight pipe except that here the direction of 5D varies 
with the axial distance z. 
From (67), we find finally 
te . ° 
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12 3 17280 4b 
Taking the Reynolds number R to be fixed and to have the value 63-3, 


(70) reduces to 


Aa: 


2 
26-38 +-2-22S ~sind. 
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This result shows that the effect of a negative S is to set up normal 
pressures which tend to keep the boundary section circular and that a 
positive S tends to make the pipe wall collapse. 

Eustice (8) has made studies of the flow of Newtonian liquids (in this 
case water) in pipes of various radii of curvature. Similar experiments 
with different non-Newtonian liquids do not appear to have been carried 
out, so that the theoretical predictions cannot at this stage be subjected 
to any observational test. 
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SUMMARY 

A new representation of the relaxation spectrum of a liquid is used in order to 
develop the theory of oscillatory flow of the most general linear elastico-viscous liquid 
in a coaxial-cylinder elastoviscometer of the type built by Oldroyd, Strawbridge, 
and Toms. It is shown that experimental results concerning dilute polymer solutions, 
hitherto interpreted in terms of a discrete relaxation spectrum, can be alternatively 
interpreted in terms of a simple continuous relaxation spectrum characterized by 
three constants. It is also shown that the ambiguities of interpretation arising from 
experiments in which the amplitude ratio alone is measured as a function of fre- 
quency could in principle be partly resolved by measuring in addition the phase-lag 
between the inner and outer cylinders of the apparatus over the same frequency 
range. 


1. Introduction 
In the first paper of the same title (1), consideration was given to the 


theory of flow of one particular elastico-viscous liquid prototype in a co- 
axial cylinder elastoviscometer of the type built by Oldroyd, Strawbridge, 


and Toms (2). The prototype considered was one which has proved very 
useful in characterizing certain dilute polymer solutions (2, 3, 4), namely, 
that represented by an equation of state of the form 


ics ‘ O\. «3 
(1 Mo) Pia 2no( I te Nee (1) 


at sufficiently small rates of strain. The notation is that of (1), ¢!}) being the 
rate of strain tensor and p’, that part of the stress tensor associated with 
change of shape of the material element. In the present paper, the earlier 
work is extended to include the most general linear elastico-viscous liquid. 
This involves a consideration of the relaxation spectrum and in order to 
deal with mobile liquid systems satisfactorily it is found convenient to 
consider a somewhat different representation of the spectrum from those 
found in the literature (see, for example, 5, 6, 7). 

Theoretical predictions are made on the basis of a number of idealized 
continuous relaxation spectra, and it is shown that experimental results 


Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 4, 1960) 
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concerning dilute polymer solutions, hitherto interpreted in terms of a 
discrete relaxation spectrum corresponding to equation (1), can be alterna- 
tively interpreted in terms of a continuous relaxation spectrum of a simple 
type requiring no more physical constants to specify it. The experimental 
method of measuring amplitude ratio alone in an elastoviscometer is in fact 
insensitive even to certain gross changes in the form of the spectrum. It is 
shown that the method can be made sensitive by measuring in addition the 
phase-lag between the cylinders. 


2. The relaxation spectrum 
In order to describe the relaxation spectrum in the most appropriate way, 
we recall first the simplest (Maxwell) model of an elastico-viscous liquid 
whose behaviour is characterized by an equation of state of the form (cf. 7) 
Pik, Pik _ 91) 
vy “eR 9 
7 
or, what is equivalent, ‘ 
2y f 


Pix ett) re'D)(t") dt’, (3) 


. 
x 


where +(= » ¢) is the relaxation time. Considering next a number of 
discrete Maxwell elements connected in parallel in the conventional spring- 
dashpot model representation, we can generalize (3) with the aid of the 
superposition principle (8), to give 

t 


ih : 

Pix > aad | et-7 Dt") dt’, (4) 
inn ¢ 

where 7; and 7; correspond to the ith Maxwell element. 

The theoretical extension to a continuous distribution of relaxation 
times has been considered by several authors (see, for example, 5, 6, 7). 
Without exception, they define a distribution function by considering the 
rigidity of the Maxwell elements as distributed over a continuous range 
of relaxation times. Since we shall be concerned with elastico-viscous 
liquids, we think it more appropriate to define a distribution function of the 
viscosity (rather than the rigidity). The distribution function of relaxation 
times (or the ‘relaxation spectrum’) N(r) is accordingly here defined such 
that \(r) dr represents the total viscosity of all the Maxwell elements with 
relaxation times between 7 and r-+-dr. The equation of state then becomes 
(generalizing (4)) 


a 


ett Dt) dt’, (5) 


3 A (7) dr 


’ 
Pa =< 
4 

. 

0 
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and if we introduce the relaxation function ’, defined by 


y N(r) | 


- 


tt’) 7 dr 
’ 


¥(t—t’) - | 
0 
(5) takes on the familiar form (cf. 7) 
: 
Pix 2 | ¥(t- t')eP(t’) dt’. 
‘ 

This equation may be considered as the most general linear equation of 
state for an elastico-viscous liquid. Furthermore, if the relaxation spec- 
trum V(r) is known for all 7, the mechanical behaviour of the material is 
determined completely, and the spectrum can be used as a perfectly general 
representation of elastico-viscous liquids (ef. 7). 

In slow steady motion (¢‘) small and independent of the time), equation 
(5) reduces to 
Je wy) | N(r) dr, (8) 

0 
and if we introduce a viscosity yo, defined as the area under the (7, .V) curve 
No | N(r) dr, (9) 
0 


equation (8) becomes Pix 2 es? (10) 


so that », can be identified with the limiting viscosity at small rates of 
shear, as observed in steady-state experiments. This is of some importance, 
since the limiting viscosity », is measured in many experiments and a 
normalizing condition (9) is immediately imposed on the relaxation 
spectrum. In the older representations of the spectrum as a distribution 
of rigidity no such simple significance can be given to the area under the 
curve which represents the spectrum graphically. 

In terms of the above representation of the relaxation spectrum, the 
spectrum associated with equation (1) is given by 


N(r) 4 











« —A,——e 


where 5 is Dirac’s 5-function. 
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3. Theoretical considerations 

We can now investigate the motion when a general liquid of the type 
detined by equation (5) is contained between coaxial cylinders and subjected 
to shearing motion. This motion is supposed brought about by imposing 
on the outer cylinder forced harmonic angular oscillations about its axis, 
the inner cylinder being constrained by a torsion wire. As in (1), all tensor 
quantities are expressed in terms of their physical components in cylin- 
drical polar coordinates (r, ¢, z), z being measured vertically upwards. Two- 
dimensional axially symmetric flow is assumed, with velocity components 
(0, rw(r, t), 0). 

The theory is very similar to that of (1), except that in this case the 
physical component of stress rd’ is given by 

: , 
ro’ = | V(t— tyes daw’. (12) 


This equation has to be solved in conjunction with the equations of motion. 
For the simple type of motion considered in this paper the important equa- 
tion of motion is that concerned with rd’, which may be written as 
eee 9,1’ ‘ 
or 2r Cw P 
yo ee (13) 
cr r ct 
where p is the density. 
Equations (12) and (13) have to be solved subject to the appropriate 
boundary conditions. A solution is possible of the form 


rp’ ref (r)e27in"t, w re F(r)e27'", (14) 
where » is the frequency of oscillation, and f(r) and F(r) incorporate appro- 
priate phase factors. Equation (12) becomes 

x 
dF ad 


W(Eje*7i"E dé, (15) 
dr 


t 
fir) ¢ 2rrinty - | W(t—t' je2rint’ dt’ — 
ar 


x 0 


and we also have, from (13), 


“(f) ~ drinpr® F. 
Eliminating F, we find 
d? d 
re _ tp 4 gt2_4)f — 0, 
( drt dr" " \s 


o..° 2 
; —2r7in —<7Ftn 
where PRs. ees. Ae 


) W(Eje-27in€ dé { (N(r) dr)/(1+-27inzr) 
0 a 
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In (1), the differential equation (17) was derived with the particular 

Se ee 2minp( 1+ 27rinA,) 

“ ccaeet i (19) 
No( 1+ 27nd.) 


This particular case can be obtained from the general expression (18) by 
substituting for V(r) from (11). 

The solution of (17) can be expressed in terms of Bessel functions of 
complex argument in the form 


f(r) aJ,(xr) + bY,(xr), (20) 


where a and / are constants to be determined from the boundary conditions. 
These are conveniently expressed in terms of the angular amplitudes of the 
two cylinders. If 6, and 6, are the amplitudes of angular oscillation of the 
inner and outer cylinders (radii r, and r,, respectively), then 

F(r,) = 2min,, F(r,) = 2in@, e*, (21) 
where ¢ is a real constant. We also have the equation of motion of the inner 
cylinder 2ar? Lf(r,) = (K—47°n?1)0,, 
where J is the moment of inertia of the inner cylinder about its axis, L 
the length of the annular gap, and A the restoring constant of the torsion 
wire. 

Since the measured quantity in the apparatus of Oldroyd, Strawbridge, 
and Toms is the ratio # of the angular amplitudes of the two cylinders, 
the solution of equation (17), subject to the boundary conditions (21) and 
(22), is most conveniently expressed in the form 

(2zpr} L) x{ J, ( xr i x5) Y,( xT} J, ( xr) | + 
-(K 47*n?- 1){ J,(ar, Yi xy) — ¥,( x )J}( wre) | (23) 

(27prir, L) x| J,( wr )¥y (ary) —Yo(ar,)J,(a7, )| 
Equation (23) is unsuitable for direct computation since there appear to 
be no tables of the Bessel functions J, and Y, for general complex arguments. 
An approximate formula, valid for sufficiently large n, is obtained by sub- 
stituting for the Bessel functions their asymptotic expansions in powers 
of (ar)~'. In this way, equation (23) reduces to the approximate formula 


l (r,\3} 3(r5—T;) 
‘ | ‘) (1 -, res 1’ s)cos(r, ra 
Uv ls | STS - 


ie | 
1) |siner -1,)a}, (24) 
x 


l5r,—3r, 15 45 
_ > XY + - 


Sr= x l6r3 x3 


~ Rye 
12875 


where 7p \("72) and 


25) 
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The above equations are of the same form as those of Oldroyd (1), but here 
x has a different significance. 

Predictions of the form of the function #(n) can now be made by con- 
sidering a number of simple idealized continuous spectra. These predictions 
may be compared with those for various discrete spectra of the type (11). 
Except where otherwise stated, the values of the apparatus constants 
used in the numerical calculations correspond to the experimental con- 
ditions of Oldroyd, Strawbridge, and Toms (2) designated by A (Cylinder 
II, Torsion wire 8). 


(a) The simple block continuous spectrum 

One of the simplest geometrical forms of the continuous relaxation 
spectrum corresponds to a rectangle in the (7, V) diagram, stretching 
from the origin. Two constants are sufficient to specify such a spectrum, the 
limiting viscosity np» (the area under the rectangle), and the length § of 
the rectangle on the time scale. We have 


' N( T) 


V(t) = 9B (O27 <8) | 
(" 





V(7) = 0 (+ > B) 
1o 


<-—_—_ 3—> 











— 7 


For this simple prototype, from (18), 
3 47*n*pB 
. Ny log(1 4 27inB) 

Fig. | contains a family of (%, n) curves for liquids of this type. The curves 
have the same general features as the published (m9, Ay, Ag) curves (1), each 
showing a peak at a frequency greater than the natural frequency n, of free 
oscillation of the inner cylinder on its torsion wire, with } passing through 
unity at a value of n almost exactly equal to no. 

It has not yet been possible to identify this new prototype with the 
behaviour of any real liquid, although it is thought that it might prove a 
useful approximation in the case of certain polymer solutions in a slightly 
higher concentration range than those considered by Oldroyd, Strawbridge, 
and Toms. 


(b) The displaced block spectrum 
An obvious generalization of the first prototype is obtained by displacing 
the rectangle along the time axis. Such a spectrum requires three constants 
5002.52 Gg 
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Fic. 1. The variation of #(n) with B when p ; 8. 


jo» By, By for its description, where the block lies between 8, and f, on the 
time scale, In this case 


7 





4 N(r) 
| 
| 
| 


| 


Br 


- a a 
pS 
N(r) 0 < B,) | 
N(r) = no/(B2—B,) (By <7 < By)). 
N(r) = 0 (r > B.) 
For general values of 8, and 8,, 
2 ; 47°n*p(B.—B;) ee 
no log {(1+-27inB,) (1+ 27inB,)} 
Figs. 2 and 3 contain two families of (3, n) curves for this new prototype. 
Fig. 2 is obtained by keeping 8, constant and allowing f, to vary, while 








(29) 
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Fig. 3 is obtained by keeping the dimensions of the block constant and 
displacing it along the time axis. The curves have, once again, some of 
the features of observed curves, but detailed comparison shows them to be 
inadequate to characterize the behaviour of the polymer solutions observed 
by Oldroyd, Strawbridge, and Toms. 


(c) The three-constant spectrum 

A more realistic prototype is one whose spectrum has a Newtonian ele- 
ment at the origin, with the non-Newtonian behaviour represented by 
a rectangular block stretching from the origin. If the Newtonian element 
has a viscosity on, (0 < o < 1), it is convenient to characterize this new 
prototype by the three constants , 0, 8, where f is the length of the 
rectangular block on the time scale. We have in this case 


N(r) 


To 





(1-0)no 


eae 








—— T 





(30) 


N(x) any 9(r) +e (0<7<p)|. 


V(r) = 0 (tr > B) 
—— AaPn®pB 
2mrinBony (1 3) No log (1 +-2minp)] 


(31) 


2 
x* 





Fig. 4 contains a representative set of (#, n) curves for three-constant 
idealized systems of this type. 

It is clear from the graphs that these prototypes lead to (3, n) curves 
which have the same features as the curves corresponding to the three- 
constant discrete spectrum (11) (1).- Since Oldroyd, Strawbridge, and 
‘Toms have found those curves useful in the interpretation of experimental 
results on dilute polymer solutions (2, 3,-4), it is of interest to attempt 
an alternative interpretation in terms of the new prototype introduced 
in the present paper. 


4. Comparison with experiment 

We confine attention to the relaxation spectrum given by (30) and 
described by the three constants 7», 0, 8. If this spectrum is to characterize 
any real liquid, n» has to be identified with the limiting viscosity at small 
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rates of shear, and o and £ can then be estimated by comparing the experi- 
mental (#, n) curve with families of theoretical curves based on the known 
p and 7, and a series of values of o and f. Figs. 5 and 6 illustrate the type of 
agreement which is possible when this new prototype is employed. It is 
clear that the experimental results, previously interpreted in terms of a 
discrete spectrum, can be equaliy well interpreted in terms of a continuous 
spectrum. It is, however, necessary to check that this interpretation of the 











1 1 


15 20 
n (c/s) 


Fic. 4. The variation of #(n) with a and B when p 1, mo 





experimental results is independent of the particular experimental arrange- 
ment chosen and we must next compare the observations on a particular 
polymer solution in a number of experimental set-ups with a corresponding 
set of theoretical curves. 

We consider for the purpose of illustration some results which have had 
a wide covering in the literature (1, 2, 9), obtained for a certain mixture 
of polymethyl methacrylate in pyridine at 25° C. This liquid, of density 
0-98 g per ml, contained 30-5 g of polymer per litre and (from experiments 
involving steady rotation of the outer cylinder) , was estimated to be 
7-9 poises. In the earlier work (2), the solution was characterized by a 
discrete spectrum with A, — 0-065 see and A, — 0-015 see. 

The values o — 0-13, 8 — 0-18 see are chosen by inspection to give the 
best fit for the experimental set-up A, and are used to construct a set of 
theoretical curves relating to five other experimental conditions. Figs. 
7-12 illustrate the agreement between theory and experiment. This agree- 
ment is thought to be somewhat closer than that obtained by considering 
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Fic. 5. The predicted relation between 3 and n in the case of an idealized 
clastico-viscous liquid for which p 0-985, ny = 3-9, o 0-25, B = 0-08; 
compared with experimental points for a solution of polymethyl metha- 
erylate in pyridine (cone. 36-4 g/l) (25° C). 








7 i 1 


40 15 20 25 
n (</s) 


Fic. 6. The predicted relation between # and n in the case of an idealized 





elastico-viscous liquid for which p= 1, m = 16, o 0-05, B = 0-79; 


compared with experimental points for a solution of polymethyl methacry- 
late in pyridine (78 per cent)/water (22 per cent) (conc. 28 g/1) (24-3° C.). 
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Fic. 7. Predicted (full line) and observed (e) relations between # and n. 
(Inner Cylinder IT, Torsion Wire 8.) 
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Fic. 8. Predicted (full line) and observed (e) relations between &# and n. 
(Inner Cylinder II, Torsion Wire 5 A.) 
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Pie. 12. Predicted (full line) and observed (e) relations between & and n. 
(Inner Cylinder IH, Torsion Wire 5 A.) 
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Fic. 15. Predicted relations between c and n. 
(Inner Cylinder II, Torsion Wire 8.) 
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Fic. 16. Predicted relations between ¢ and n. 
(Inner Cylinder III, Torsion Wire 8.) 
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the (n», Ay, A.) model of Oldroyd (1), but the improvement is certainly not 
sufficient to allow one to accept the new interpretation to the exclusion 
of the old. Our conclusion must be that, as far as the experimental results 
of Oldroyd, Strawbridge, and Toms are concerned, the simple continuous 
spectrum defined by the three constants (np, o, 8) is adequate to explain 
the behaviour of dilute polymer solutions in a simple shearing motion at 
small rates of strain. This means that Oldroyd’s original interpretation in 
terms of a three-constant model is not unique. 

The findings of this section lead one to the conclusion that the experi- 
mental method of measuring amplitude ratio alone is insensitive even to 
certain gross changes in the form of the spectrum, and we are led to con 
sider whether the method can be made sensitive by measuring the phase- 
lage between the cylinders as well. 


5. The study of phase-lag 

The phase-lag between the two cylinders was not measured in the experi- 

ments of Oldroyd, Strawbridge, and Toms, although it is in principle 

possible to measure phases and amplitudes simultaneously; and later 

designs of the same instrument do in fact include facilities for measuring 

both amplitude ratio and phase-lag (10, 11). 

No new theory is required, since the phase-lag of the inner behind the 

outer cylinder is ¢ in equation (21), given by 
r 7pr? L) x] oJ,( wr Pi xy) Y,( ur edi ( u's) | ! 
(A 4rPn?— 1) A (ar Vy (org) V(r A (ore) |} (32) 
((2aprir, L) ald,(ar,)¥{(ar,)— ¥2(ar,)4,(07,) |} 


Making the same approximations as in the derivation of (24), this becomes 


f ri \o{ 3(7.—T7,) 
F (1 : ! sjeos(r, ry)a 
v fa | syn = 


"0 


1 


87's oy l6r3 


| Lor, — 3r 15 


3 $5 . | a 
(. : 4g] [SIM(T2 r;)x,. (33) 
x? Sr= x 12875 Ny . | 
) 


0 
KMquation (33) is not valid for small values of n, and in order to determine 
the limiting value of ¢ at small frequencies, we need to use the series expan 
sions for the Bessel functions in equation (32). For very small values 
of n, equation (32) becomes 


(34) 


which implies that —$ — 2 Sb ; (35) 


This result is independent of the relaxation spectrum and is true for all 
linear elastico-viscous liquids. 
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Figs. 13 and 14 contain some representative theoretical (c, n) curves. 
Fig. 13 relates to discrete spectra of the type (11), corresponding to the par- 
ticular (3, m) curves contained in Fig. 1 of (1). Fig. 14 relates to simple 
block continuous spectra of the type (26), corresponding to the particular 
(3, x) curves shown in Fig. 1 of the present paper. It will be noticed that ¢ 
passes through zero at frequencies almost exactly equal to the natural 
frequency mo, a result which is independent of the type of relaxation spec- 
trum considered and may be compared with the approximate coincidence 
of ® — Land n — np in the amplitude-ratio curves, 

Figs. 15 and 16 contain, for two sets of experimental conditions, the 
theoretical (¢, n) curves for two spectra, indistinguishable in the amplitude- 
ratio experiments of Oldroyd, Strawbridge, and Toms. The difference 
between the curves at high frequencies (as much as 20°) would certainly 
be detectable in phase-lag experiments, so that measurements of phase-lag 
could easily resolve the ambiguities arising from amplitude-ratio experi- 
ments. We conclude, therefore, that such measurements ought to be made, 
together with amplitude ratio, if a complete characterization of elastico- 


viscous liquids is required. 
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UNIFORMLY STRETCHED PLATES SUBJECTED 
TO CONCENTRATED TRANSVERSE FORCES 


By ARNOLD D. KERR 


(Institute of Mathematical Sciences, New York University) 
{Reeeived 21 August 1959. tevise received 7 January 1960] 


SUMMARY 

This paper contains a study of the behaviour of isotropic clastic plates of various 
hapes subjected to uniform tension in the plane of the plate and loaded transversely 
by concentrated forces. 

The deflexion w of the plate is governed by the partial differential equation 
\Aw—(N D)Au 0, where N is the tension intensity per unit length and D is 
the flexural rigidity of the plate. The fundamental deflexion function (Green’s 
function for an unbounded domain) is determined and used in connexion with the 
method of images to construct solutions for plates of various shapes, simply sup- 
ported along their boundaries. 

Solutions are obtained for (a) the wedge-shaped plate with opening angle a = 7m 
m 1, 2,3,...) and (6) the reetangular plate. It is shown that the rectangular corner 
plate, the infinite and semi-infinite strip, can be obtained as special cases. The 


rectangular corner plate is discussed in more detail, 


1. The mathematical formulation of the problem 

THe formulation is based on the assumption of the so-called small 
deflexion theory of plates. For the assumptions and derivations the reader 
is referred to the standard texts on the subject, e.g. (1). The deflexion 


normal to the surface of the plate w due to a laterally distributed load q 
is governed by 


AAw— x2Aw 7 (1) 


where Ais the Laplace operator in two independent variables and «? — ND. 
Here N is the tension intensity per unit length and D — FA3/12(1—,?*) 
is the flexural rigidity of the plate. 

The moments, shear forces, and the reaction distributions (along the 
simply supported edges parallel to the y- or x-axis) are formed as higher 
derivatives of w as follows: 

O2y ae aw 
M, 4 M, = Di! 
2 v n 2 
cy” 
Sen 
—D(i—p)- 
excy 


Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 4, 1960) 
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ox CXCY™ 
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2. The fundamental solution 
The fundamental solution of (1) is defined here as a solution axially 
symmetrical in nature due to the action of a concentrated force P (q — 0). 


S ld ‘i 
alk K"1w 0. 
dr? eh 


w— Cylnr+€, +0, 1(«r) +, Kg(«r), (6) 


For this case (1) reduces to 


ld 
dr? = a) 


Its general solution is 


where /,(«r) and A,(«r) are modified Bessel functions of order zero of first 
or second kind respectively. 

Assuming w == Oatr — Oand due to the condition at r -> « the deflexion 
surface reduces to the form 


uw Ollnr + Ko(«r))]. (7) 


To determine C we consider the vertical equilibrium of a centrally 
loaded element separated from the rest of the plate by a cylindrical surface 
of radius r and obtain p 

C= = 0 (8) 
2aN 
Hence, the fundamental deflexion of an elastic plate subjected in its plane 
to hydrostatic tension and laterally to a concentrated force P is 
P . 
w= poy line + Kg(«r)|, (‘)) 
where r is the distance from P to the investigated point of the plate. 
Through expansion of Ag(«r) it can be shown that (9) contains the character- 
istic singularity (1, 87/))r?Inr of a plate in bending due to a concentrated 
force P. 

Each term on the right-hand side of (9) separately satisfies (5), since 
Inr is a solution of Aw — 0, and A,(«r) is a solution of (A—«*)w = 0, Thus 
considering (9) it can be concluded that the deflexion of a uniformly 





464 A. D. KERR 


stressed plate subjected laterally to P consists of the deflexion of a uni- 
formly stressed membrane subjected to P, w = —(P/27N)Inr, and the 
solution of (A—«®)w — 0 which allows for the ‘flexural rigidity’ of the 
plate. It is in fact the In-term of the flexural rigidity expression which 


makes w finite at r 0. 


3. Solution of boundary value problems 

For plates bounded by simply supported straight boundaries solutions 
can be obtained using the ‘method of images’ (1). A more general method 
suitable for plates of arbitrary shape and arbitrary degree of fixity is based 
on the fact that the deflexion can be represented in the form w = wy,-+-2w,, 
where w, is the fundamental deflexion and w, are functions which satisfy 
A\Aw, — 2? Aw, 0, adjusting the w, values at the boundaries to the pre- 
scribed ones, The exactness of w depends then mainly on how accurately w, 
can correct w, at the boundary to the prescribed values. It should be noted 
that the method of images is a special case of the above. 

In the following, plates of various shapes will be treated by the method 


of images. 


4. Wedge-shaped plates. The rectangular corner plate 

Solutions for wedge-shaped plates can be obtained in closed form by 
arranging an even number of forces P on a circle about the tip of a plate. 
For a plate with an opening angle « 7m (m 1, 2,...) subjected to a 
force P in an arbitrary point inside the plate 


2m 


- ¥ (—1)"*Ylnr, + Ay(xr)], (10) 
2aN n=1 f 


where the 7,,’s are the distances of the forces P to the point of the plate 
under investigation. For the reetangular corner plate m —= 2 (Fig. 1). 
(10) reduces to 

"13 


+ Ky(«r,) — Ag(«r) + Ag(«rs)— Ag(Kr,) |, (11) 


P lin 
7N | 


a) ry 
where 
1 = Vt(Z—ApP+(y—Yo)*}, Ta = VU(T— Tp)? +(Y+ Ho), < 

rs = V{(r+29)?+(y+y,)"k, ry = V{(x +29)? + (y—Yo)*}. 

It can be shown that the boundary conditions w = 0 and Aw = 0 along 
the supported edges :.ve satisfied. 

To find the influence of the ratio «2 ND upon the deflexions, (11) is 
reduced for the case when P acts at a point on the bisector radius (at r = ry) 
to an equation of the deflexion curve along that line. The numerical evalua- 


tion of this equation for different values «? is shown in Fig. 2. 
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Fic. 1. Force system for the rectangular corner plate. 
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Fic. 2. Corner plate. Defiexion curves along the bisector radius for 
different values of x. 1, 5. 


It is interesting to note that whereas for the case NV + 0 the deflexions 
at infinity tend to zero, for the particular case NV 


0 the deflexions at 
infinity are infinite. 


Because of the assumptions made in deriving (1) a negative concentrated 
corner reaction appears in the analytical solution in addition to the reactions 
distributed along the supports: 


R = 2D(1—p) =a ie (13) 


Cxey y=, 
Hh 
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Substituting (11) into (13) and setting x = 0 and y = 0 we obtain 


R—tM=H)Pl_ 220% = Xoo 


K,{ x+y). 14 
7 \«?(az 4 y2)? (a2 +y2) atv (Zot Ho) _— 


J 
This is the corner reaction for an arbitrary position of P inside the plate. 
The numerical evaluation of R as function of x? for the case x) = yy = a 
is shown in Fig. 3. 
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‘orner plate. Corner reaction R as function of x?. 
be 0-3, a 10, 
The distributed reactions along the supports are obtained substituting 
(11) into (4). For example the reactions at the support along the y-axis are 


Vv) P(2(1 ef) Lol 3(y—Yo)?- xo], xf 3(Y+Yo)?—2x5 
K? | x2-+-(y—yp)?]8 [234 -(y4 Yo)? | 


pw) Ky(«y {x3 +-(y—Yo)?})—1 J+ 


3 (Lp) Ka(xy (25 + (¥+-Yo)?}) —1]4 


Yo)” 
ro(¥ Yo)" K. {2 2) 
“(x2 L(y—y aL 3(Ky {20+ (Y¥—Yo)*})— 

0 ‘ 
LAY > Yo)” 
(x24 (y+Yo 2 |! 





Kyle (23+ (y-4 wo)|| | 


To illustrate the influence of the stretching forces upon the reaction distribu- 
tions, graphs of [V) |.) for « = 0-3 and different values of «x? are shown 
in Fig. 4, for the case x, = y, = a. We find that the reactions for various 
values of x? are bracketed between two limiting cases, namely, that of 
x* — 0 (plate in bending) and the one corresponding to «x? = oo (stretched 
membrane). It is to be noted that the maximum value of [V,],_, does not 
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occur at the point on the boundary closest to P, not even in the case of 
stretched membrane. It can also be observed that the areas enclosed by 
the reaction distribution for different x? values (but for the same value 
of the force P) are not equal. This is due to the fact that the corner reaction 
R varies with changing x* and because of the vertical equilibrium, for 
decreasing R the area enclosed by V has to decrease since 


x 


R+P—2 | [¥,],-0dy = 0. (16) 











Fig. 4. Corner plate. Reaction distribution along y = 0 or x = 0. 
K xo or 0, ---- «? l, ' jo» st -—*? ibs» ed 0-3, a 10. 


5. The rectangular plate 

A system of forces as shown in Fig. 5 periodically repeating in the positive 
and negative directions of the 2- and y-axis produces a series of deflexion 
surfaces each of them equivalent to that of a rectangular plate simply 
supported along its boundaries and subjected at an arbitrary point (x9, 4) 
to a force P. 

The deflexion surface is 


> — tes , 
— I I TnmPnm 
w=-s y ni — ; _ 
ons m=— Oo n= —a@ | TamPnm 


aie Ky(xr, nd + Ko(«T am) i Ky(xp,, a) a K,( Pan) ’ ( l 7) 


where , ae 
y{(a@—2na— xy)? + (y—2mb—y,)*} 


y{(w7— 2na— xy)? +(y— 2mb+ yo)} 
(18) 
v{(w—2na+-x9)?+(y—2mb—yo)*} 


- 4 {(w@—2na+-xy)? + (y—2mb-+ yo)?} 





7) Force 


The bending moments are derived by substituting (17) into the first two 


relations of (2). 
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system for the rectangular plate 
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2 4 
<a \ K : Tram 


2na)?—(y—Yo—2mb)* (x 4 ty 2na)*?— (y+ Yo — 2mb)? 
nm Pnm 


2na)*— (y+ Yq— 2mb)* (1 -}- pe) 
Kk |r 
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K,(«r 


nm 


1 n m) 
p nm 


, , _ , , 
Ky (KT an) Ky (KPa) tT , K (Pn) T 


’ nin Pn m nm 


| u(r -ty — 2na)*-+-(Y— Yo 2mb)*| 5 ler 


nm ) 
nm 


7 


| pa(ar tty — 2na)*+ (Y— Yo 2mb)?) i ( ’ 
» 2 K nna 


Tam 
| u(x a“ suey" (Y} Yo eT) Kin.) ; 
Pnm 


| ye(r-} ary 2na)? + (y+ yy 2mb)*] )- ro 
Pom 7 





, 
| 
nm ) | 


Kk KT ny) 1 Ky( Kvn) } Ky lk py, m) Ko KPa) 


(21) 

The sufficient condition for validity of term-by-term differentiation of 
an infinite series in two variables is satisfied since (17) as well as its first, 
second, and third derivatives with respect to 2 and y are continuous fune- 
tions of cand yinthe whole domain (except the second and third derivatives 
ats y — Owhere they have a singularity). 

For solutions expressed in terms of infinite series, the problem of con- 
vergence, as well as the rate of convergence, is of considerable importance. 
From the theory of intinites series it is known (2) that if the terms of a 
ad Mg ses (22 
be of alternate sign, it is necessary and sufficient for the convergence of 
the series that 

u, ° U,., for every value of n, and lim u, — 0. (23) 

n+ 
Because } P — Othe solution obtained using the method of images forms 
(for each point on the plate) an alternating series. It can be shown 
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that the condition in (23) is satisfied and that the resulting series are con- 
vergent. The validity of condition (23) for the cases treated here can also 


be proved by considering the physical picture. Namely, if a stretched 
plate is subjected to a concentrated force P, the effects (like w, M, V) 


are localized and will decrease with increasing distance from P. Astheterms 
of increasing n represent the effects of forces P acting at increasing distances 
from the point under consideration, these terms will decrease with increasing 


15 401 
ie | U 


0 





i ia 


0:08 | 
p Var]ac=0 
P 
Fie. 6. Infinite strip. Reaction distributions along x = 0 or x = a for pp = 0°25. 


2 s 1 1 
we? xr or 0, ~ eal, «? os» — K? too @ 10, 





n, and this proves physically that condition (23) is satisfied for the cases 
under consideration. 

In representing the deflexion surface due to P by means of a Fourier 
series, as done in (1), 310, we deal with a force system similar to the one 
represented in Fig. 5. In the case of the sine series, with n increasing, each 
additional term improves the approximate representation of the ‘con- 
centrated’ forces, all of which are already acting on the infinite plate. 
The double sine series satisfies automatically the condition of simply sup- 
ported boundaries. In the case of a representation of the deflexion surface 
by means of fundamental deflexions, (17), each additional term represents 
an additional ‘complete’ force P applied on the infinite plate, which con- 
tributes towards the formation of simply supported boundaries. 

From what has been said above, one can anticipate the very slow con- 
vergence of the formulae for the bending moments, expressed as double 
sine series, in the neighbourhood of P, since it is rather ‘inconvenient’ for a 
sine or cosine series to represent these strong irregularities. There will be a 
quicker convergence for moments derived by means of the fundamental 
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solutions, since in the vicinity of (x9, y,) the main part of the moments is 
represented by the fundamental solution of P at (x9, yo), and the contribu- 
tion of moments caused by the other forces P is diminishing rapidly with 
increasing distance from (29, Yo). 

It is to be noted that for m = 0 and n = 0 the expressions obtained for 
the rectangular plate reduce to that of a rectangular corner plate. he 
summation of the terms p, 9 and 7, .forn — 1, 2,... only, yields the solution 
for a semi-infinite strip. Summing the r, ,’s for n — 1, 2,... only, we obtain 
the solution for an infinite strip. In Fig. 6 the distribution of reactions is 
shown for the infinite strip of width a, subjected to a force P at (0, 4a). 
Also, here it can be seen that the reactions for various values of «? are 
enclosed between the two limiting cases: that of a stretched membrane and 
the one of a plate in bending. 
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SUMMARY 
In this paper the Green's funetion is derived for a layered fluid of finite depth 
ubject to a veneral boundary condition at the bottom. The limiting form of the 
Green's function as the depth is increased indefinitely is considered, and a physical 
interpretation of the results given 
Phe relation between the limiting form of our problem and the same problem set 


vith a radiation condition at *% is examined \ particular example is worked in 
deta 


The form of the solution due to a transient pulse is noted; also some aspects of 


the problem of murmerical solution A yveneral formula for the group velocity ts 
derived 


1. Introduction 

lit basie problem to be discussed here is the solution of the wave equation 
ina medium in which the velocity of propagation ¢ depends on the depth 
coordinate z alone. In these circumstances the wave equation can be 
written 


(1) 


or, in the slightly more general form, 


( fle P| ep 
p) 


a) he ' O2Rp oz} fc(z)}* at? 


which allows for changes in density p with depth. The latter equation is 
amenable to the same treatment as (1) though we shall have little cause 
to use it. An exception occurs when p is allowed to have step discon 
tinuities corresponding to layers of immiscible thuid of different densities 
as in section 5, 

We consider first solutions which have a simple harmonic dependence 


on time. This leads us to consider the equation (Helmholtz equation) 


veep i 0. (la) 


(0(2)5* 


Together with this equation we consider the boundary conditions: 


(a) radially outgoing wave (or radiation condition) as r—— . (474 y?) > & 


Quart. Journ, Mech. and Applied Math., Vol. XIII, Pt. 4, 1960) 





PROPAGATION OF SOUND IN A LAYERED FLUID MEDIUM — 473 
For our purposes it will be satisfactory to say that f(v) satisfies a radiation 
condition at if 
S(v) + Ae +g(r), and g(v) +0, v > x, 


This condition is a particular case (with g(v) suitably restricted) of the 
Sommerfeld condition f 


lim wv! = Lif 0 


vs at 
which has been employed as a sufficient condition in proving uniqueness 
theorems for infinite regions (see for example, Stoker (1)); 
(b) Py, pressure release surface; and 


CP. 


cz 


(c) cosaPe yt sina 0. 

Condition (¢) contains P(/1) (pressure release bottom), and 
PH) ez 0 (rigid bottom) as special cases so that it is physically 
appropriate. To understand its mathematical significance we introduce 
eylindrical coordinates in equation (la) and then separate the variables. 

sv putting P — O(r)d(z) we obtain the pair of equations 


V2 04 k70 0 


and a ' (, al ke \4 (), 
( 


dz {e(z)}? 


where 4? is the separation constant. 

Thus we see that (¢) is the most general condition at 2 — // which 
depends solely on the values of 4 and dd dz at z H, and which makes 
equation (3) self-adjoint on (0,//). Thus equation (3) has a discrete 
sequence of real eigenvalues k2 + — & and eigenfunctions ¢,, as its possible 
solutions. Boundary condition (a) implies that we must choose the Bessel 
function //)*) as the solution of (2). 

The choice of k k* as the separation constant leads to a more con- 
ventional form of the eigenvalue problem [equation (3)], and to the modi- 
fied Bessel function of the second kind K (kr) as the solution of (2) again 
perhaps the more conventional choice. We have adopted our notation 
to conform with that used by Pekeris in his paper (2) in order to be able 
to make direct comparisons with his work. 


2. Construction of the Green’s function 


To find the Green’s function for (1a), that is the solution due to a point 
source, we assume that it can be expressed in terms of the Green's 
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functions of (2) and (3); this suggests that its form is 


$,,(2)bn (Zp) 
fhe)” 


where we have to determine f(k2), and where we have written 


G(r: Z,Zp) oe hi S HY (k,, r) 
— 


n=0 


r = .[(x—azp)*+(y—yp)*)]. 
We will assume in this section that the ¢,, are normalized so that 
H 


| ¢. dz = 1. 
6 
Now G(r: z,zp) has the defining property that, if 
Vb + [w? fe(z) 2 |® — R(F)Z(z), 
then (P) — 20 | | G(r: z,zp)R(F)Z(z)F dFdz, 
where P denotes the point with coordinates (xp, yp, zp) and 7 = , (4? +-y?). 
Substituting for G(r: z,z,) we get 


@(P) = ¥ — Ani | He? (k, RUF ar | Pnl2)On(2p) 72) G2\ 
@ \J 


f(k) | 


n 

> @.(P). 
ike n\ 

n 

Now we have 


V2 ® (P) k2@, (P)+ Rirp) | Pol2)Pol2r) Z(z) dz 
: J fk) 
since — }iH\?'(k,,r) is the Green’s function of (2), and 
oD, (P) aw? 


ez (C(Zp)§* 


- db (P) ke @®,(P), 


i 


since ¢,,(zp) satisfies (3). That is, 
TDP) +P) = Rirp) [ Zi) Y Pn@bn GN ae 
(f°) fe(zp)}2 (f°) re) | we f(ke : 
R(rp)Z(zp). 
5 $,,(2)4,,(Zp) 
cone {(k) 


which implies that (4?) 1. Thus we have 


Hence 3(z—Zp), 


G(r: 2, zp) bi SY HP (ky, T)bn(2)bn(Zp)- (4) 
0 


n 
This series converges for all values of its argument which lie within the 
volume considered provided that we interpret v(—k?) as —1\k'. The Bessel 
function Hi? then decays exponentially for r > 0 as k2 + —o. 
We here introduce the term ‘mode’ to describe the individual terms in (4). 
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Those modes with k2 < 0 are sometimes called ‘evanescent’ to distinguish 
them from the propagated modes (k2 > 0). 

The argument we have used above can be made rigorous for the sequence 
of integrals /,, defined by 


* k 
I, = | Z(z)\ > Pul2)ba(20)| dz. 


It defines the generalized function 6(z—-z,p) provided Z(z) is suitably 
restricted (see, for instance, Lighthill (3)). 


3. Spectral analysis of equation (3) 
In this section we consider the limiting form of the Green’s function 
given in equation (4) as the depth is increased without limit (H -» x). 
We assume that ¢(z) - positive limit < 0. Therefore w?/{¢(z)}* — limit 
a® » 0. For convenience we write equation (3) as 


d *h ; { w 
f 
‘ 


+ 


dz* 


a?+a? Kl 0, 


2 
or 4 {q(z)+-Ajd 0. (3*) 


2 
c(z)}? 


Provided that q(z)\dz < a, then (3*) hasa fundamental set of solutions 
M(z) > et (z -> w), 
and N(z) > e-**4* (z > 2) 


(these being the leading terms of asymptotic expansions valid in a neigh- 
bourhood of z =). In particular N(z) satisfies a radiation condition 
at «. If |q(z)| is not integrable on (0,00) then we find asymptotic expan- 
sions which have a more complicated form but which behave similarly 
for large values of z. To simplify our arguments we will assume that |q(z)! 
is integrable. 

The characteristic equation for equation (3) (the equation whose zeros 
as a function of k* give the eigenvalues) is 

M(0) N(0) 

dM(H) cosal(ll)+-slae dN(H) = 0. 


cos aM (H)-+-sin « dz dz 


We now show that, if —q,,,, < A < 0, then A(H) takes discrete values 
which tend to definite limits independent of « as H -- «. This follows 
because 

(a) for finite H the problem is self-adjoint with, therefore, a discrete 
spectrum, and hence only a finite number of eigenvalues in this range; 
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(b) M and N are integral functions of A for finite (fixed) values of 2. 
see (4), chapters | and 7: 


(c) we have M(z)>ca (z—> 0), 
and N(z)>0 (z—-0), 


for A in the stated range. Hence the characteristic equation for large H 
can be written 
M(0)[cos aN (#1) + sin a}d N(H) dz} | 


N(O) : 7 
cos x~M(H)+sinajd. M(H) dz} 


and the magnitude of the right-hand side can be made arbitrarily small 
by taking H large enough. Thus, as H > x, the roots of this equation 
tend to the roots of V(0) 0: and, as (0) is an integral function of A, 
the roots of NV(0)— 0 cannot have a finite limit point. By a similar 
argument we can show that the eigenfunctions also tend to definite limits. 
and that the difference between any eigenfunction and its limit function 
can be made arbitrarily small on (0, 7) by taking H large enough. 

In limit 47 + x values of A > 0 lie in the continuous spectrum of (3). 
To show this we have to prove that the number of eigenvalues between 
any two fixed values of A increases without limit as H/—> x. 

That this must be so follows from the oscillation properties of the eigen 
functions the nth eigenfunction being characterized by having exactly 
n zeros in the interior of (0,77). From the asymptotic forms for WM and 
we see that these oscillation properties must be similar to those of sin vAz 
for A 0 and = large enough. The change in the number of zeros of 
sin <Az on (0,77) within a given range of A tends to « w H->+ x, hence 
so must the number of eigenvalues of our problem which lie in this range. 

We can give a physical interpretation to these results if we associate 
with each value of A a direction 6 in accordance with the rule 

sin (2). 
C(z) 

This allows us to associate a real ray with each mude for k*? > 0; and 
these rays have the properties required of them by geometrical optics. 
For example we have 


w C(2,) 


‘ WwW ° 
sin 6(2z) sin 4(z,), 
C(2,) (2) €(2,) 


, of). 
with sin @(2,) D sin O(2). 


Thus the refraction law of geometrical optics holds with @ interpreted 
as the angle between the ray and the axis of z. 
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When —dmax < A < 0 there is, for each ray, a critical depth at which 
it is refracted back towards the surface. At depths greater than this the 
corresponding mode decays exponentially. These are modes that are 
formed by the constructive interference of signals describing refraction 
paths in the fluid. They are thus independent of the nature of the bottom 
provided it is deep enough. 

Those modes with A > 0 (a? > k* > 0) correspond to rays which suffer 
reflections at the bottom. We can thus interpret these modes as being 
formed by the constructive interference of bottom reflected signals. When 
H — x& we see that there is no interference mechanism, and all signal 
paths are equally to be preferred. To this extent there is a similarity 
between the continuous spectrum and an outgoing wave (see also the 
example of section 5 and equations (6) and (7) of section 6). 

The form of equation (4) is not valid in the limit // © if equation (3) 
has a continuous spectrum. We therefore modify (4) by introducing the 
function o(f*) which is constant except at the eigenvalues k? — k2 when 

H H 
it jumps by 1/ | d2 dz (that is, o(k2-+-0)—o(k2—0) = 1 ‘fg dz). We can 
£ 


0 0 
then write the Green’s function as the Stieltjes integral 


k§=~« 


G(r: z, Zp) Hg Hi? (kr)b(z, k)db(zp, k) da. 
ke 


\s in section 2 the convention , (— k*) -ik must be observed. 

The limit of o(4*) as H —» & is again a step function for /* in the discrete 
spectrum. However, o is continuous (even differentiable) for * in the 
continuous spectrum. For proofs of these results see, for example, 
Coddington and Levinson ((4), chapter 9), or Titchmarsh ((5), chapters 
2-3). 


4. Pekeris’s construction for the Green’s function 

An alternative method of obtaining the result (4) is to construct 
G(r: 2,2») as the sum of the residues of a contour integral in the complex 
k-plane. These poles are at the zeros of the period equation as defined by 
Pekeris in (2). We now derive the period equations in the sense of 
Pekeris 

(a) for the problem as we have stated it; and 

(6) for an outgoing wave condition at z 20. 

Following Pekeris we begin by attempting to find a solution which has 
a point source singularity at (r = 0, <= zp). Thus we seek a solution 
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which behaves like 


e e—thz-zp % 
Jy(kr) 0 -kdk as ,/{r?+(z—zp)*} > 9, 


Q? = w?/{e(zp)}?#—k*. 
We try the solution P = | Jjkr)d(z,zp)k dk. 
0 
As ¢ satisfies the second-order differential equation (3), we see that the 
conditions 


(1) f(z = zp—O) = O(z = zp+9) 


and 
i See eee ¥.s 


dz dz 
are sufficient to specify the solution uniquely. 
We put ¢,(z) = A, M(z)+-B,N(z) for z < zp, 
and $,(z) = A,M(z)+B,N(z) forz >2z 


ope 
The boundary conditions then give 


d, (2) AS N(0)M(z)—M(0)N(z)} 


and 
$,(2) 


B!| NUA)cos x— } dM (HH) 
\ 


sin « |.M(z) | MU eos —~ sin x N(2)! 


dN(H) 
2 dz { 


The matching conditions then give 
(1) A{N(O)M(z,)— MO)N(zp)} 


B | N(H)cos «+ = (#) sin »| Mer) 
| M(H )cos a+ aM(H) 


and 


(2) Aly) M(0) 


ade ( 


dN(zp)| . 
dz | 

_IN(A) dM (zp) 
dz ¥ dz 


OS x 


=? 


adM(H) . dN(zp)) 
sin |" 4 


| M(H)ecos « 


These equations can be written in matrix form as 


Vx = b. 
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The period equation is then 
det(V) = 0. 
Explicitly we have 


M(zp) N(zp) N(0) N(H)cos «+ — sin « 


Y A] [0 
dM (zp) dN(zp) dM(H) la - Bi 


dz dz —M(0) —M(H)cosa— -— q 


~ 


We see that the first matrix is the matrix of the Wronskian of the 
differential equation (3) and hence is non-singular as M and N are a 
fundamental set of solutions; and that the second matrix is the matrix 
of the determinant whose zeros (as functions of k?) are the eigenvalues. 

Thus we can always solve for A and B unless ‘? is an eigenvalue of 
equation (3). This establishes the equivalence of the zeros of the period 
equation and the eigenvalues in this case. 

In the case of an outgoing wave at z = 2% we have 


$,(z) = B, N(z), 
and the matching conditions lead to the matrix equation 


M(zp) N(zp) 


A N(0) = O|[A 0 
amen) aXer|| “wey a|Le] = [of 
dz dz ‘ 


We see that the characteristic equation is 
N(0) = 0, 


which is the limit of the characteristic equation for the boundary condi- 
tions (a) as H -» « provided that —q,,,, < A < 0. Thus we again recover 
the modes associated with the refracted ray paths. The outgoing wave 
condition is necessarily set on an infinite range so that it is of interest to 
ask what, in this case, corresponds to the continuous spectrum. Pekeris, 
who set the outgoing wave condition in his work, found, in addition to 
the refracted modes, an integral term to which he gave the name of 
‘branch line integral’. In the next section we reproduce as an example 
one of the problems solved by Pekeris, and show that, in this case, the 
branch line integral is identical with the term to which the continuous 
spectrum gives rise. 

Our approach in this section has been somewhat more general than 
that of Pekeris’s original work and is due in part to Haskell, see Ewing (6) 
this reference also contains a full summary of Pekeris’s results. The 
use of contour integration to derive eigenfunction expansions is described 
in detail in Titchmarsh (5). 
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5. Ona problem solved by Pekeris 
When the density p is allowed to vary with z, equation (3) becomes 
é/! dd\ Vf w? 
dz 4 * plfe(z)}? 
This equation is again self-adjoint; and our previous arguments can be 
applied to it. However, in this case, the orthogonality relation becomes 
ul 


7“ 
b; >; d> 
pP 
and equation (4) becomes 
. p z)d Zp 
G(r whi H'2(ke, ry Pn! Pni=y) 
$ <~ p(=p) 


i ) 


Also, if p is discontinuous at certain depths, then we must have the 
additional condition 

ldd 

continuous 

pe 
satisfied at these px are the Weierstrass. Erdmann corner 
conditions: see Coura bHhiibert (7). They correspond to continuity of 
pressure and aormatl ve acity at the layer interface, 


We will consider tiv 


continuities at I 


pecial model where both p and ¢ have step dis 


jare otherwise constant. Our boundary conditions 


are P 0 at Oana // 
for the case // f 


This problem was solved by Pekeris (2) 
ind then used for interpreting the results of explosive 
sound transmissions in shallow water over a fluid bottom 

The eigenfunctions are given by 


42,2 in layer (1), 
and | ke) Bsingt(77 2) im layer (2) 
N \C9 


The eondition that the matching conditions be satisfied leads to the 
characteristic equation 


ind, D singd,(/7 D) 


| 
QQ, cosQ, D MO cos Q,(I 1D) 
2) Pe 


For w* ¢ i wo c2 the speetrum is discrete in the limit // 


The characteristic equation in this case is 


tanQ, D ' tanh (2, (// 1D) 
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Py 442 


This equation has discrete solutions if it has solutions at all, The limiting 
| 4 


which tends to tanQ, D as H+. 


form of the eigenfunction in layer (2) is given by 
lim p,sinh Q, (ff 2) pa e'teiD-s) 


H +x ‘bt, cosh Q, (H—D) i2, 


9 


The spectrum is continuous in the limit H —- if w/e > k*. We can 


write the characteristic equation as 
tanQ, 7 tanQ, D py Sd, 


tand2, D. 
1-+tanQ, //tanQ, D pry, — 


We see that, if 14-tanQ, /tanQ, ) — 0, then it vanishes again for 


O, O,4AQ,, AQ, 1+ A ya) 


Now exactly one eigenvalue lies in the range (Q,,Q,), Thus the number 
of eigenvalues lying in any finite range (Q,,Q, 4 AQ,) tends to infinity with 
AQ, Ho, that is with J. 

It is of interest to verify that the eigenfunctions are orthogonal on 
(0,1) with respect to the weight function 1p (as they should be). This 
orthogonality property eluded Pekeris, a fact on which he comments at 
some length. 

Qur eigenfunctions have the form 

p,sin@2, z 

(2, cos@2, D 

p,sind2,(/7 2) 

(2, cosd2,(/1 dD) 
» " 

We put db; dz, and I, . db; p; dz, 
Pi J Pe 


0 Db 


db 


in layer (1), 


and d in layer (2). 


Then 
1 p 
j 1 
VOB OF, 0,42), cos, D costa, D” 


“1Q,,cosQ,, DsinQ,, DQ, sinQ,, DeosQ, D; 


(22, 003, 02,,02,, cos, (HI " Dyoowtd,,(H D) ’ 

“1Q,,cosQ,,(H— D)sinQ, (1D) Q,, sinQ,,(H— D)cos Q,,(H D)}. 
4(D) -4(D) — 4(D)-&(D) 
QF, —Q3; k? — k} 
Ii 


I, 


’ 
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d,(D)—4,(D) $;(D)—4d,(D) 
and . ; ee all 
O38, —- 23; ki—kF 


Hence i+h=90 (t #9). 
We see that the role of the weight function is essential to the argument. 
In this case the spectral function can be constructed; and thus we can 
give explicitly the limiting form of equation (4). However, we cannot 


use the representation 

p, sin Q,( Hf —z) 

0, cosQ,(H D) 
for the solution in layer (2) as this expression is indeterminate in the limit 
H x. A suitable form is 


PQ, cosQ, DsinQ,(2—D) 


For large values of /7 we have 


|_O( D/A). 


\lso, for large values of 77, the distance between adjacent eigenvalues 


so} > by 
is given b (AQ) H a+O(1 A). 


Thus we have 


-O(D #H),. 


\o(Q,) HT H/p2Q3 sin2?Q, D+ p3.Q2 cos?Q, D 
AQ, 7 2 ; 025 cos?Q, D 


dalQ ‘ ()20)2 anger 
ia lo {2,) 9 827 O5 cos*Q, D 


AQ, 7 p02 sin2Q, D+ p3Q? cos?Q, D 
Let denote the limiting form of that part of (4) which belongs to the 
continuous spectrum, Assuming that zp belongs to layer 1, we have in 
this laver 
128in 22, zp 
5 QF cos*Q, D 


k dk. 


For terms in the discrete spectrum we have 
1 | p; sin*Q, z 2 | de. | l (p3 7 
Py \Q2F cos *0), D| py\ 02, |? 


202, D 2 | 
dz 
| 
D 
(D sin2Q, D| Po 
QF cos*Q, D\2 402, a 2 02, 3 
92 3 


Pi lO, D—sinQ, DcosQ, D+ P27, 
# 1 cos*Q2, D\ Py “0, 
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Since tanQ, D= — pahdy 
Pp; 82, 
we get 
I ye dz - pI —— 
p 2023 cos?Q, D 
«{Q, D—sinQ, DcosQ, D—(p, py)? sin’?Q, DtanQ, D}. 
Thus the limiting form of (4) is 
G(r: 2, 2p) 
hi < H2(k, ) | Q, sinQ, z —_ am a if 
od Q, D—sinQ, DcosQ, D—(p, py)*sin®?Q, DtanQ, D 
mice H2\ ker)? p2 92, sinQ, zsinQ, 2), 


— —— k dk. 
p7 $25 sin?Q, D+ p§QF cos? Q, D 


The form of the solution given by Pekeris for this problem is 
4aGi(r: 2, Zp). 

We can explain the factor 47 if we note that we have effectively adopted 
the form (1 47R)e-‘” (R? — r?-4-(z—z,))*) for the singularity at z = zp 
while that used by Pekeris is there equivalent to (1 Rje-t@”*, 
6. Propagation of a transient pulse 

For completeness we sketch how our solutions can be modified to 
describe the propagation of a transient pulse. Let q(w) be the Fourier 
transform of this pulse. Then the required solution is given by 


x 
. 


I ; 
P(r,z,t) = = giw)G(r: z, 2p, wel dw, 


. 
x 


where we consider our Green’s function as a function of the para- 
meter w. 

Let us assume that our observations are made at a sufficiently great 
distance ry so that we may use the asymptotic representations of the 


Hankel function ie 
H® k iii = ~4 imtkhyr) 
eee de a 


x 
. 


> ix | ae ; _ _ 
sii S7r p | no) | ;}Pn(=)bn(=0) ila+k, r—wl} yy, 


We then have 


In practice the major contributions to these integrals are near the points 
of stationary phase of the exponential terms (see, for instance, Lamb (8) 
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where the ¢, are smoothly varying functions of w). In this case we can 
make use of the asymptotic formula 


x 


( d(xjel™ dx 


x 


v7 b(z) 
v( df"(%) ) 


etl f(z)+ 17} 


where we take +}7 according to the sign of f”(%), and # is the (here 
assumed unique) point of stationary phase (that is f’(7) = 0). 
In our problem the points of stationary phase are given by 


1 
~-fot—k.#) = 0, 


n 


Iw 
° dw 
that is by r t Cf. 
; . dk " 


n 
where "(', is the group velocity of the nth mode. 


The equation pa OCS 


determines (for a fixed value of r) the arrival times of the different 
frequency components of the nth mode (here considered as a function of 
w, and characterized by the property that ¢, vanishes n times in the 
interior of (0,7) for all values of w). Thus we see that the arrivals at a 
given point at a given instant are made up of components from different 
modes which have the same group velocity. In general each of the arrivals 
will depend on distinct values of w (generally different for each mode) 
unless (', as a function of w, is stationary for this value of w. When this 
occurs one would expect some degree of reinforcement. This has been 
called the ‘Airy phase’ by Pekeris. The method of stationary phase 
breaks down for the Airy phase. Pekeris has given a formula which is 
here applicable. 


When d{"C_} dw + 0, the solution is given by 


; “ ] dine ’ 1 . 
P ~ ler ya go) a 1 by (2)by (2p) t{hn 


where the summation is taken over those values of w such that 


r=", ana & > @. 


The computation of P from formula (5) requires, in general, the 
numerical solution of equation (3) for, at any rate, the first few eigenvalues 
and eigenfunctions for a set of values of w within the desired frequency 
range—any other values can then be filled in by interpolation. This is 
a perfectly practical task using modern computing methods. To compute 


di"(\ dw it is first necessary to know "C), and this can be determined 





PROPAGATION OF SOUND IN A LAYERED FLUID MEDIUM 


numerically by the formula H 
| bi(z)de 
nC, = 8, yg —__—__. 
| [ew {e(2)}*]bi(2) dz 
provided we know &? and 4¢,,. : 


To derive this we differentiate equation (3) with respect to w. This 


© feb) (a pallbe _fog dba 20 | 
d2\éw)” \tez)}® "few "dw — fe(2)}"I 


gives 


pn: 


As the boundary conditions hold for all w we must have the right-hand 
side orthogonal to ¢,, over (0, H) for a solution to be possible. Hence 


1e(2)}") 


H 
| a .-* een ne, 
F 


From this expression the formula for the group velocity follows imme- 
diately. When & is in the continuous spectrum, we have 


C,=k/—® 
. Haan 


H 


(as | ¢2 dz is divergent in limit H = o for this value of k). 
0 
Putting k = {w c(z)\sin @(z) we see that 


C22 = ¢(x)sin 6(x). (7) 


That is, the group velocity ultimately corresponds to the horizontal 
component of the velocity along the ray corresponding to the given value 
of k. This provides further evidence for the identification of the continuous 
spectrum with an outgoing wave. 

Formula (6) is perhaps most valuable if difference methods are used to 
solve the eigenvalue problem (3) for & in the discrete spectrum. See, for 
instance, Fox (9), or Kopal (10). Some recent developments of the 
methods described in (9) have been given by Osborne (11). 

I am indebted to a referee for bringing to my notice some recent work 
on geometric diffraction theory. He has suggested that these methods 
would be of value in obtaining numerical results. My impression is that 
a direct construction of the geometric approximation is useful only before 
the shape of the pulse becomes distorted by dispersion (that is at short 
ranges), while (5) is most valuable after dispersion has separated out the 
various mode terms. Thus, the two approaches would seem to be comple- 
mentary. In one of the references (Friedlander (12), chapter 3, sections 
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8-9) some results linking geometrical acoustics with the theory of harmonic 


wave trains are given. Although this theory serves to unify the two 


approaches, in practice (see (12), chapter 6) this often means that a 


programme similar to the one we have adopted is used. 
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SUMMARY 

Finite-difference approximations to second-order differential equations are 
generally less satisfactory when terms involving first derivatives are present. We 
here discuss some methods of approximation which are aimed at improving the 
accuracy in such cases, They arise out of a closer investigation of a method originally 
proposed by D.N. de G. Allen and R. V. Southwell for dealing with a partial 
differential equation which occurs in the field of viscous fluid motion. 

The paper falls into three fairly distinet parts. In the first we examine the basic 
approximation and we are then led to consider whether equally satisfactory results 
cannot be obtained in certain cases when first-derivative terms are not present. 
Both ordinary and partial differential equations are considered up to this point; 
we then nuke an extension of the basic theory of more especial interest in the case 


of ordinary ditferential equations, 


1. Introduction 

IN a recent paper, Allen and Southwell (1) have described a new treatment 
by relaxation methods of a partial differential equation which occurs in 
the tield of viscous fluid motion. With slight changes in notation the equa- 
tion is 


0, (1) 


| oC res) eC ops cL 


R\ cx? cycr) cy oxy 

where Cis the vorticity, % the stream function, R the Reynolds number, 
and « and y orthogonal coordinates. In the numerical treatment of (1) the 
function yy may be regarded as a known function of 2 and y (it is in fact 
found from a second equation to be solved jointly with (1)). For con- 
venience we write cb cor —- f(a,y) and ob cy — g(x,y). Then the standard 
relaxation treatment, using a square mesh of side / with points numbered 
as in Fig. L, is to replace the ¢-derivatives at 0 by their finite-difference 
equivalents (Southwell (2)) which leads, with neglect of terms of order h‘, 
to the approximation 


(5 go h)or (1+ bfyh)l, 1 (’ 


y y » 1\, 
LR T Soh |b (1 Ay hl, (1 i \eo v. 


Rk, 
(2) 
The aim of the relaxation process is to satisfy (2) at each nodal point 0 to 
some degree of accuracy acceptable as error (residual error). This can be 
made as small as we please for a given mesh size but the solution possesses 


[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 4, 1960) 
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also the more inherent error due to the finite-difference approximation. If 
f(x,y) and g(x,y) are numerically large or vary rapidly this error can be 
serious except on a very fine mesh. Allen and Southwell suggest that this 
is so in their case and have devised the following alternative approximation 
to (1), aimed at overcoming the difficulty. In the neighbourhood of 0 the 








44 


functions f(x,y) and g(x,y) are replaced by their nodal values fy and go 
and (1) is separated into the two equations 


¥ 


Rg” RA | 


= fo . é | 
y~ cy 


where A is constant. The solution of the first of (3) is 

dol — Ax + P+-Qekoor, (4) 
where P and @ are integration constants. Using (4) as an expansion along 
the x-line 301 we can dispose of P and Q and express A in terms of nodal 
values of C, thus: 


oa v 


Joi (Cy Co) Rach _; S3 Cos Ah(e~ Raer — 1), (5) 


Dealing similarly with the second equation of (3), which likewise yields an 
expansion to be used along the y-line 402, we obtain 


fot (Ss Coe Joh Ce Ce Ah(et* 1). (6) 





ON FINITE DIFFERENCES 
Equating the values of A in (5) and (6) gives the equation 
Jo h{(l, — ole a" + L,— Co} /(1 —e- Far) — 

Soht(G- Cy )ef” { Fi Co} (1 —eSeh) 0, (7) 
which Allen and Southwell find to be a superior nodal approximation to (2). 
It is clear that any superiority of (7) over (2) must lie in the approximation 
it makes to the first-derivative terms in (1), for it is easily shown that the 
two formulae are identical when these terms are absent. It is primarily this 
question that we investigate in this paper, i.e. we attempt to discover the 
reason for the superiority and to point out any advantages in this type of 
approach which is undoubtedly new in finite-difference theory. 


2. It is convenient to start with the ordinary differential equation 
C+ pla)’ +-q(a)o = r(x), (8) 
where primes denote differentiation with respect to x, and subsequently to 
make the extension to partial equations. There are, in any case, some 
points of special interest for ordinary equations. The foregoing process 
would replace (8) along the line 301 by the equation 
o" + pol +o = To (9) 


To1. tper( P cosh kyx+-Qsinh ky x), (10) 


do 


where ko = Gp2—q)!. (11) 


so that 


Eliminating P and Y in terms of nodal values of [ gives 


» 

y , 2r, 

L, clvoh 4 Ly etek 20, cosh ky h : "(cosh bpyh—coshkyh) — 0, (12) 
0 


which is sufficient to define a set of nodal approximations to the one- 
dimensional solution of (8). Equation (12) is valid whether ky is real or 
imaginary. It holds also when q, — 0 (although (10) must now be modified) 
and if we take the appropriate limit of the last term we obtain 


" 2rah . 
C, etre +L, eto — 22 cosh dpgh 0" sinh Jpyh — 0, (13) 
Po 

which is the same as (5) when applied to the first equation of (3) (where A 
would be given if this were an independent ordinary equation). Now the 
standard approximation to (8) derives from the central-difference formulae, 
in usual notation (Fox (3)) 

Ure (5 — 189+ gud> — 1 4Hd7-+-...)Cy (14) 
and hy = (8? — 84+ J8°— 459+... )Lo, (15) 
and the customary three-point formula 


(L+-4poh)l,+(1—4pyh)l,—(2—h?qy)lo—h?ry = 0 (16) 
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is obtained by neglecting differences beyond the second in (14) and (15). 
The reason that it is sometimes unsatisfactory when the first-derivative 
term is important is evident from the much slower rate of convergence of 
the series (14) as compared with (15) (which is usually highly satisfactory). 
\ possible way of getting over the difficulty is to remove the first-derivative 


term in (8) by the usual substitution 
C Ye-ifp de, (17) 
where the new dependent variable satisfies the equation 


(1S) 
ind oe a Lp ° (19) 
There are sometimes objections to this approach. For example, the new 
variable may not be well behaved over the range of integration: or if p(x) 
undergoes large variations, the right hand side of (1S) may vary too 
rapidly for satisfactory numerical treatment. But both these objections 
can be removed if we apply (17) only locally over the interval 301 and 
moreover onli porarily, 1.e. having applied it we return at once to the 
variable €. Moreover, when we do this we tind that the resulting approxima- 
tion bears a very close relationship to Allen and Southwell’s, particularly 


when it is applied to (1) which is, rather fortuitously, a very special case. 


3. The basic approximation 
We suppose the points 3, 0,1 (Fig. 1) to be respectively 2,— 4, x9, and 


r, - hand let (17) be detined in the interval (2, —A, 15+ 4). We may choose 


0 
any value of the indetinite integral so for convenience we take the value 


which vanishes at 2 ry. 1.e. let 


0 


u(r) L | pdr (20) 


so that the corresponding form of (1S) is 
Y"+¢4¥ = re’. (21) 
Neglecting terms in /' the standard three-point approximation to (21) is 
y, +-¥,—(2—h*ho)¥,—h2ry = 0 (22) 
and eliminating pivotal values of Y using (17) we have finally 
Ga {,e"?—(2—h*ho)Co— hr (), (23) 
If pis given as an analytical function of x, and u can be found exactly, then 


(23) is the final basic approximation, Otherwise we must expand p about 


ry in the usual way, Le. 


P | °° 
p(r) Pot (X ~LoQ)Pot 5; (% Le)"Pot-+-5 (24) 
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and a first approximation to u would be 
u(x) = bpo(r7—Xy), 
giving the corresponding form of (23) 
C, elvoh + Lye bred —(2—h2bg)lg—h*rg = 0 (26) 


We can express /*¢, in (26) in the usual finite-difference form, i.e. 


. ‘ — - 
h*dy h*(qy Po) AtD P3)- (27 


At this stage we can usefully compare (26) with (12); the exponential terms 
ineach compare exactly and if we expand the hyperbolic functions in (12) in 
powers of / as far as terms in /* we get, approximately, 


C,¢ Ho C3 ¢ bpoh _() h*qy Lh? pely hr, (), (28) 


Equations (26) and (28) agree only if p(x) is invariant with respect to x, 
ie, the last bracket in (27) is identically zero, But there ean be little doubt 
that if » varies even only moderately with 2 this term will be unportant ; 
its magnitude could easily exceed that of the other terms in /7d, (e.g. in the 
numerical example which follows). However, in obtaining (26) from (23) 
we have neglected a similar term of magnitude $/?), in approximating to 
vu, and u, and if, as we suggest, the last term in (27) should not be ignored, 
then we should also retain this term, i.e. we include the term in pg in (24) to 
give (according to the trapezoidal formula) 

h h ™ 

| Pot Py); us 4'Pe + py). (=) 
More accurate formulae can be obtained by taking into account more terms 
of (24) (e.g. the three-point formulae obtained by assuming parabolic 
variation of p over the points 3, 0, and 1). 


Exram ple | 
A very simple example will suffice to illustrate the point. The function 
' 1 —erfr, (30) 
where erf is the error function (2 vz) | ¢ “ dx, satisfies the equation 
0 


ver 


¢°+22f’ = 0, (31) 
with boundary conditions ((0) 1, U(x) © It is recorded in Table 1 
atintervalsh — 0-4 and correct to four decimals; also recorded are approxi- 
mate ¢-values satisfying respectively the standard approximation (16) 
(with py 22. Jo — To — 0), equation (13), and equation (23) (with u- 
values calculated from (29), which gives exact values of the integrals in 
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this case). The approximate solutions were obtained by the usual relaxa- 
tion methods, i.e. by liquidating residuals, defined by the left-hand side 
of the appropriate equation, using relaxation patterns. The results ade- 
quately confirm the points we have made so far, which we may summarize 
thus. Although the order of the error term in (16) and (23) appears to be 
the same, the latter is more satisfactory because it eliminates the effect 
of the slowly convergent difference series associated with the first-derivative 
term. Equation (12) (or (13) if appropriate, as in this example) has the same 


TABLE | 





Eqn. (13) Eqn. (23) Eqn. (30) 





y y 


~ 4 S 





1°0000 1*0000 
0°5726 o°5716 
0°2583 0°2579 
00595 0°03597 
0°0238 0°0237 
0°0047 0°0047 
0°0007 0°0007 
o*000! O°0001 














effect but omits an important term arising from the variation of p(z). 
We can show that (12) is certainly correct to the same order of accuracy 
as (16). Expanding the exponential and hyperbolic functions in (12) we 


obtain 
(i- bpoh Cy ~(1 spoh)ls (2 h?qo)Lo : h*r, T bpeh(g, i {.— 20) T 
i Lpanr(t, l3) O(h*) (32) 


and the last two terms on the left-hand side of (32) involve 82, and pdf, 
respectively and so, by (14) and (15), are each of order h*; thus (12) and 
(16) ditfer by terms of order h* (but whose magnitude may be large if p, is 
large). 


4. To deal with the partial differential equation (1) we proceed in a 
similar manner to that used by Allen and Southwell, separating it into the 


two equations y ‘ 

G2 

> _ Ra(x y) — RA(x,y) | 
cr cr 


(33) 


and x,y) A(x, y) 


cy cy 


Along the x-line 301 we have g(x, y) = g(x, yo) and A(x,y) = A(x, yp), 
where 0 is (29, ¥)), and the solution for Z is a function of x alone. Along the 
y-line 402, f(x,y) = f(xy, y), A(x, y) = A(x, y) and (isa function of y alone. 
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Applying the process of the last section to each of (33) in turn and elimi- 
nating Ay (= A(x», Y»)) from the two resulting equations we obtain as an 
approximation to (1) at 0 the equation 


a 
Ro 


» 
em +f, evs) +{, ets + byers—(2 +l) to = @, (34) 


where x 
u(x) = —3R | g(x, yo) dx | 


To 


v(y) = 4 | f(x,y) dy | 


Yo 
and P, (2) - ; (5 — }(f2+ Rg). (36) 
- 0 


To obtain some comparison between (34) and (7) we can proceed as in the 
previous section. Equation (7) can be written in the equivalent form ob- 
tained by combining two equations of type (13) (one in the x-direction and 
one in the y-direction) and we can expand each of these in the form (28) 
(with appropriate changes in symbols). Thus combining these last two 


we obtain 

fo = 9 
(37) 

as a close approximation to (7), i.e. with neglect of terms in /*. To find 

expressions for pivotal values of u and v we proceed as in the last section 


by expanding g(2, ¥) along the 2-line 301 and likewise expanding f(2x,, y) 
along the y-line 402. Substituting the corresponding first approximations 


S Fa al 
4 Rgoh | : $Rgoh) 3) ® bfoh | »—4foh __ ls an ei 2) R. 2 
Cg ¢ ) f,¢ lye \ k 4 (f go); 


found in this way into (34) this equation becomes 
» 
iKaoh +f, eh Raoh) 1 7, ehfoh. tye tin—(2 rs phe) 0, (38) 


which we can now compare with (37). We find that the equations agree 
exactly by virtue of the fact that f(a, y) = Gb x, g(x, y) = &b/cy, so that 
the first two terms on the right-hand side of (36) cancel each other identi- 
cally. From the point of view of the numerical analysis this is rather fortui- 
tous; in the general case f and g may make an important contribution to 
(36) which our analysis suggests should not be ignored. The other general 
point which emerges is that even in Allen and Southwell’s case a sub- 
stantial improvement in accuracy should result from using the trapezoidal 
(or an even more accurate) formula to estimate pivotal values of u and vr; 
this would involve an almost negligible increase in labour. 
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5. This point leads us to consider whether Allen and Southwell’s process 
could itself be modified to allow of this better approximation. As it stands 


the method replaces (8) by the equation (9); but if we first apply the trans- 
formation (17) to (8) and then approximate to the resulting equation (21) 
by the equation y’+¢,Y ry (39) 
we would obtain 
es 
y+ ¥,—2Y, cos dj, h — ~"°(1—cos di h) (40) 


Po 


as an approximation to (21) at each nodal point and hence 


: 2r : 
C,€"*— 20, cos b3 h *(l—eos di h) = 0 (41) 
Po 


as a nodal approximation to (8), « having the significance of (20). This 
formula is quite different from either (12) or (23) although it agrees with 
the latter as far as terms in /*. It has two clear advantages over (12) in that 


we can extend ul and “ 


:' , bevond the first approximations u, — }poh, 


U. \p,h (to which (12) is of necessity restricted) and, moreover, it 


T ikes prope 1 


per account of the term involving p’ (which occurs ind). Applying 
of the equations (33) in turn and combining the results 


uusly used, weobtain asa nodalapproximation to (1) the 


osh vi h) 


cosh 
hey? | 





Smear a 





Pe | ee | ee a ee ae are 


ON FINITE DIFFERENCES 


6. Some applications to a class of eigenvalue problems 
We shall consider the one-dimensional problem governed by the equation 


£"+{Ap(a)—o(x)}Z = 0 (44) 


with extensions to its two-dimensional counterpart 


9 \ 


V°0+-{Ap(x, y)—o(x, y)}f = 0 (ve eal 2 (45) 
€ 


cy? f 


where in each case p > 0 throughout the domain of the problem and A is a 
parameter possessing a discrete set of eigenvalues, with a corresponding 
set of modes for ¢, when { is required to satisfy specified boundary con- 
ditions (e.g. € vanishes on a simple closed boundary in the case of (45) or 
at the end points in the case of (44)). Both equations have physical impor- 
tance in many fields. 
For (44) the usual difference approximation at a nodal point 0 is 
C1 +f3—(2—h7qg)Lo 0 Ap—a) (46) 


and to find an approximation to a given mode we have to solve (46) for a set 
of nodal ¢-values together with the corresponding eigenvalue, If an initial 
approximation to the mode of (46) is known, relaxation methods can usually 
be employed satisfactorily for this purpose, the eigenvalue A being estimated 
at any stage by the Rayleigh quotient which corresponds to (46), i.e. 


a Poso 
where the summations extend over all (-values, including the end-points. 
This method is particularly satisfactory when ¢ Hat the end points and 
the mode considered has no mternal zeros corresponding tothe smallest A 

of (46) depend upon A and converge to the true solutions 
reased mndetinitely. Ina reeent paper Bolten and Seo 
tytatedd the convergence m the specnml case , ! ia) 


thee ne chmenssonal Sehrodmger equatron wat! 


m thats Vib) chetpte thee apeper 





also the more inherent error due to the finite-difference dpproximation. If 
f(x,y) and g(x,y) are numerically large or vary rapidly this error can be 
serious except on a very fine mesh. Allen and Southwell suggest that this 
is so in their case and have devised the following alternatWe approximation 
to (1), aimed at overcoming the difficulty. In the neighbourhood of 0 the 


9s 








4 


funetions f(7,y) and g(v,y) are replaced by their nodal values fo and go 
and (1) is separated into the two equations 


ie RA | 
cr cr : 
- ; (3) 
ak. fo 2 A 
cy cy 
where A is constant. The solution of the first of (3) is 
Jo$ = Ax+ P+ Qekoor, (4) 


where P and Q are integration constants. Using (4) as an expansion along 
the z-line 301 we can dispose of P and Q and express A in terms of nodal 
values of f, thus: 


Jot (S; —Lo)e~ Rae" 4 C3—Co} = Ah(e- Ror), (5) 


Dealing similarly with the second equation of (3), which likewise yields an 
expansion to be used along the y-line 402, we obtain 


Sot(Ss ae Coetor = ty Se Co : Ah(ef* i 1). (6) 











Equating the values of A in (5) and (6) gives the equation 

Jo h{(l,—Lo)e #9" +L, —L,} (1—e-#4)— 

Soh{(S, —Cy)efor T C.— Co} (1 —eSeh) 0, (7) 

which Allen and Southwell find to be a superior nodal approximation to (2). 
[t is clear that any superiority of (7) over (2) must lie in the approximation 
it makes to the first-derivative terms in (1), for it is easily shown that the 

vo formulac are identical when these terms are absent. It is primarily thi 
question that we investigate in this paper, i.e. we attempt to discover the 
reason for the superiority and to point out any advantages in this type of 
approach which is undoubtedly new in finite-difference theory. 


2. It is convenient to start with the ordinary differential equation 
C4 play +qiaryt © r(x), (8) 
where primes denote differentiation with respect to 7, and subsequently to 
make the extension to partial equations, There are, in any case, some 


points of special interest for ordinary equations. The foregoing process 
would replace (8) along the line 301 by the equation 


C+ pol + doh — To (9) 

so that C46 t( Posh kor +Qsinh ky x), (10) 
Ye 

where kyo ~ (Lp? —qo)?. (11) 


Eliminating ? and Q in terms of nodal values of ¢ gives 
» 
. 2r, 
C, elves Lye tee — 20, cosh kgh — ~—* (cosh }pyh—coshkyh) = 9, (12) 
qo 


which is sufficient to define a set of nodal approximations to the one- 
dimensional solution of (8). Equation (12) is valid whether k, is real or 
imaginary. It holds also when q, = 0 (although (10) must now be modified) 
and if we take the appropriate limit of the last term we obtain 


» 
C, etveh +f, etre — 20, cosh 4 pyh — mA sinh {py yh = 0, (13) 
Po 


which is the same as (5) when applied to the first equation of (3) (where A 
would be given if this were an independent ordinary equation). Now the 
standard approximation to (8) derives from the central-difference formulae, 
in usual notation (Fox (3)) 


hLy = (wS—Aud?+ gyud>— houd?+...)Ly (14) 

and h2ly = (8?@— 484+ f5°— 4959+...) Lo, (15) 
and the customary three-point formula 

(1+4poh)f, + (1—dpoh)l3—(2—h*qo)Lo—h?rg = 0 (16) 
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is obtained by neglecting differences beyond the second in (14) and (15). 
The reason that it is sometimes unsatisfactory when the first-derivative 
term is important is evident from the much slower rate of convergence of 
the series (14) as compared with (15) (which is usually highly satisfactory). 
A possible way of getting over the difficulty is to remove the first-derivative 
term in (8) by the usual substitution 


C= Ye-tpar, (17) 
where the new dependent variable satisfies the equation 


Y"idY — ress az, (18) 


, 


and d(x) = q—hp'—}p?. (19) 
There are sometimes objections to this approach. For example, the new 
variable may not be well behaved over the range of integration; or if p(x) 
undergoes large variations, the right-hand side of (18) may vary too 
rapidly for satisfactory numerical treatment. But both these objections 
can be removed if we apply (17) only locally over the interval 301 and 
moreover only temporarily, i.e. having applied it we return at once to the 
variable €. Moreover, when we do this we tind that the resulting approxima- 
tion bears a very close relationship to Allen and Southwell’s, particularly 
when it is applied to (1) which is, rather fortuitously, a very special case. 


3. The basic approximation 

We suppose the points 3, 0, 1 (Fig. 1) to be respectively r5—h, x), and 
ry 4 and let (17) be detined in the interval (2)—h, 2,+h). We may choose 
any value of the indetinite integral so for convenience we take the value 
which vanishes at x4 — 2p, i.e. let 


a 
r 


u(r) \ ( pdx (20) 


r 
To 


so that the corresponding form of (18) is 

Y"+oY = re". (21) 
Neglecting terms in /* the standard three-point approximation to (21) is 
Y,+Y¥3—(2—h?¢,)¥,—h?ry = 0 (22) 

and eliminating pivotal values of Y using (17) we have finally 
¢,e"'4 C3e" (2 h*bo)Lo— hr, 0. (23) 
If p is given as an analytical function of x, and u can be found exactly, then 
(23) is the final basic approximation. Otherwise we must expand p about 


x — 2X, in the usual way, i.e. 


y ] — 
p(x) Pot+(X%—2p) Po 4+ = (e- Lq)*Pot--:s (24) 
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and a first approximation to u would be 
u(x) = 4po(x—2p), (25) 
giving the corresponding form of (23) 
L, elven +f, e-tvoh (2 _h24,.)L,—h2r, = 0. (26) 


We can express /*¢, in (26) in the usual finite-difference form, i.e. 
‘ ‘ s, | 7 
h*dy = h?(qo—14 Po) — | (Pa Ps): (27) 


At this stage we can usefully compare (26) with (12); the exponential terms 
in each compare exactly and if we expand the hyperbolic functions in (12) in 
powers of / as far as terms in /4 we get, approximately, 


¥y 


C ebro + Ly @ beet —(2—h?qg + ph2pp)lo—h*ry = 0. (28) 
Equations (26) and (28) agree only if p(x) is invariant with respect to x, 
i.e, the last bracket in (27) is identically zero. But there can be little doubt 
that if » varies even only moderately with x this term will be important; 
its magnitude could easily exceed that of the other terms in h*d, (e.g. in the 
numerical example which follows). However, in obtaining (26) from (23) 
we have neglected a similar term of magnitude }/?), in approximating to 
u, and u, and if, as we suggest, the last term in (27) should not be ignored, 
then we should also retain this term, i.e. we include the term in pg in (24) to 
give (according to the trapezoidal formula) 


h h 
, Po + Pr)s Us - 4 (2 + Py). (29) 
More accurate formulae can be obtained by taking into account more terms 


of (24) (e.g. the three-point formulae obtained by assuming parabolic 
variation of p over the points 3, 0, and 1). 


Example | 


A very simple example will suffice to illustrate the point. The function 


4 1—erfr, (30) 
where erfx is the error function (2 yz) | e-** dx, satisfies the equation 
0 
t’+ 220’ = 0, (31) 


with boundary conditions {(0) = 1, ((%) = 0. It is recorded in Table 1 
at intervals h — 0-4 and correct to four decimals; also recorded are approxi- 
mate ¢-values satisfying respectively the standard approximation (16) 
(with po — 2x9, Yo = To = 0), equation (13), and equation (23) (with w- 
values calculated from (29), which gives exact values of the integrals in 
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this case). The approximate solutions were obtained by the usual relaxa- 
tion methods, i.e. by liquidating residuals, defined by the left-hand side 
of the appropriate equation, using relaxation patterns. The results ade- 
quately confirm the points we have made so far, which we may summarize 
thus. Although the order of the error term in (16) and (23) appears to be 
the same, the latter is more satisfactory because it eliminates the effect 
of the slowly convergent difference series associated with the first-derivative 
term. Equation (12) (or (13) if appropriate, as in this example) has the same 


TABLE 1 





Ieqn. (16) Iegn. (13) Eqn. (23) Eqn. (30) 
c 4 c 4 


"0000 1°0000 1*°0000 1°0000 








0°5663 0°5726 o°5716 
o°2514 O72583 0°2579 
o°0554 005495 0°03897 
o°0215 0°02 35 0°0237 


0°0042 0°0047 0°0047 





0*0000 o000f 0°0007 0°0007 


0000 o*Oo0o0! o"ooo! o"ooo! 














effect but omits an important term arising from the variation of p(x). 
We can show that (12) is certainly correct to the same order of accuracy 
as (16). Expanding the exponential and hyperbolic functions in (12) we 
obtain 
(1+ dp hl, + (1 —dpgh)ly—(2—h?qg)lo—h?rot bpp h*(C, + Cg— 205) 4 

i AL pahr(l, l3) O(h*) (32) 
and the last two terms on the left-hand side of (32) involve 6%, and pdf, 
respectively and so, by (14) and (15), are each of order h*; thus (12) and 
(16) ditfer by terms of order A* (but whose magnitude may be large if p, is 
large). 


4. To deal with the partial differential equation (1) we proceed in a 
similar manner to that used by Allen and Southwell, separating it into the 
two equations ' ay 

RA(x,y) | 
(33) 


and > -+- f(x,y) — A(x, y) 
q y 


Along the x-line 301 we have g(a, y) = g(x, yg) and A(a,y) = A(x, yo), 
where 0 is (9, yp), and the solution for is a function of 2 alone. Along the 
y-line 402, f(x,y) = f(x, y), A(t, y) = A(x, y) and Z is a function of y alone. 
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Applying the process of the last section to each of (33) in turn and elimi- 
nating Ay (— A(x», yy)) from the two resulting equations we obtain as an 
approximation to (1) at 0 the equation 


» 
- 


] 
pire +-{,e")-+4 fe" } Cen (2 + R 19) 0, (34) 


oa 


where - 
u(x) —4R Y(X, Yo) dx 


to 


7] : 
v(y) = 4 | fro) dy | 
Vo 

ond oe =) (Z _1( 2-4 Ryi). (36) 

2\er]y 2\ey/o 

To obtain some comparison between (34) and (7) we can proceed as in the 
previous section. Equation (7) can be written in the equivalent form ob- 
tained by combining two equations of type (13) (one in the x-direction and 
one in the y-direction) and we can expand each of these in the form (28) 
(with appropriate changes in symbols). Thus combining these last two 
we obtain 


fo = 0 
(37) 


9 2 
i Rgoh | £,, ¢4 Haoh) | Cy ehfoh | Cye bfoh __ 24 at ; (f2 | Ry?) 


as a close approximation to (7), i.e. with neglect of terms in h‘. To find 
expressions for pivotal values of u and v we proceed as in the last section 
by expanding g(4, ¥9) along the x-line 301 and likewise expanding f(2», y) 
along the y-line 402, Substituting the corresponding first approximations 
found in this way into (34) this equation becomes 


R 


(Ae LD ee ae a (24 7 


» 

ph Pe)bo = 0 (88) 
which we can now compare with (37). We find that the equations agree 
exactly by virtue of the fact that f(x, y) — O&p/ex, g(x, y) — eyb/ey, so that 
the first two terms on the right-hand side of (36) cancel each other identi- 
cally. From the point of view of the numerical analysis this is rather fortui- 
tous; in the general case f and g may make an important contribution to 
(36) which our analysis suggests should not be ignored. The other general 
point which emerges is that even in Allen and Southwell’s case a sub- 
stantial improvement in accuracy should result from using the trapezoidal 
(or an even more accurate) formula to estimate pivotal values of u and rv; 
this would involve an almost negligible increase in labour. 
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5. This point leads us to consider whether Allen and Southwell’s process 
could itself be modified to allow of this better approximation. As it stands 
the method replaces (8) by the equation (9); but if we first apply the trans- 
formation (17) to (8) and then approximate to the resulting equation (21) 
by the equation 


y’+4,¥ =f, (39) 


we would obtain 


y; | Y; 2¥, cos d}h at cos di} h) 0 (40) 


0 


as an approximation to (21) at each nodal point and hence 


2f,cos d) h “ro -cosdih) = 0 (41) 
Po 
as a nodal approximation to (8), w having the significance of (20). This 
formula is quite different from either (12) or (23) although it agrees with 
the latter as far as terms in h*. It has two clear advantages over (12) in that 
we can extend ~, and uw, beyond the first approximations u, = $poh, 
Us ipyh (to which (12) is of necessity restricted) and, moreover, it 
takes proper account of the term involving p’ (which occurs ind). Applying 
this process to each of the equations (33) in turn and combining the results 
in the manner previously used, we obtain asa nodal approximation to (1) the 
equation 
fe" —2¢, cosh af h) R(1—cosh ag h)-4 


By(C. 6+, e"'— 2¢, cosh Bi h)/(1—cosh Bih) = 0, (42) 


where pK Lai | 
H/o 


and By (e). + Hfo J 


Here wu and + have the significance of (35). 


(43) 


[t would appear that (42) is a more correct application to (1) of Allen and 
Southwell’s type of process than (7) and the amount of initial computation 
in using it would be substantially the same. But in deciding its merit 
relative to (34) an essentially new point is involved. It is whether we may 
in general expect the approximation (40) to represent (21) more accurately 


than customary finite differences. There is some evidence that this may be 
so for if é and re“ were invariant with respect to x, (40) would represent (21) 
exactly. The new problem is quite a distinct one since the first-derivative 
term is now absent. In the second part of this paper we investigate this 
question, with interesting results, in relation to an important class of 
second-order equations. 
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6. Some applications to a class of eigenvalue problems 
We shall consider the one-dimensional problem governed by the equation 


{+ {Aplar)—o()}Z = 0 (44) 


with extensions to its two-dimensional counterpart 


-~ os 
TL Mptey—olnmye=0  (V= +5) (48) 
ore =6cy” 
where in each case p > 0 throughout the domain of the problem and A is a 
parameter possessing a discrete set of eigenvalues, with a corresponding 
set of modes for f, when { is required to satisfy specified boundary con- 
ditions (e.g. € vanishes on a simple closed boundary in the case of (45) or 
at the end tal in the case of (44)). Both equations have physical impor- 
tance in many fields. 
For (44) the usual difference approximation at a nodal point 0 is 
0.+¢,—(2—h?qy)ly = 0 (q = Ap—o) (46) 
and to find an approximation to a given mode we have to solve (46) for a set 
of nodal ¢-values together with the corresponding eigenvalue. If an initial 
approximation to the mode of (46) is known, relaxation methods can usually 
be employed satisfactorily for this purpose, the eigenvalue A being estimated 
at any stage by the Rayleigh quotient which corresponds to (46), i.e. 


Nhe 3 > ((o4 L3—2lq)Lo— -h? “OoS Co} 
> po 


(47) 


where the summations extend over all ¢-values, including the end-points. 
This method is particularly satisfactory when ¢ = 0 at the end-points and 
the mode considered has no internal zeros (corresponding to the smallest A). 
The solutions of (46) depend upon A and converge to the true solutions 
of (44) as / is decreased indefinitely. In a recent paper, Bolton and Scoins 
(4) have investigated the convergence in the special case p(x) = 1; (44) 
then becomes the one-dimensional Schrédinger equation with o as the 
potential function. They show that if A — A(h) denotes the approximation 
obtained from (46) to a given true A of (44), then 

Ath) = A+-vh?+ O(h) (48) 
and that if ¢ = 0 at the end-points and o is bounded in (0, 1) and non- 
singular at the end-points, then v < 0, i.e. for small enough h, A(h) < A. 
We shall show that this result is valid for the more general case of (44) but 
first let us consider the modifications introduced if we apply Allen and 


Southwell’s process to (44). The nodal approximation is 


0,+f,—2f,cosgih = 0 (49) 
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which we can write approximately as 
Ci4+-Cg—(2—h?qo4 hhtqi)oo 0, (50) 


for small enough A, so that the order of difference between (46) and (49) 
is at once apparent. We can obtain a Rayleigh quotient corresponding to 
either (49) or (50); one can be obtained from the latter by assuming the 
term h4q? C,/12 to be temporarily independent of changes in A. This gives 


Pais y ¢; (24 Whrgilotlo hog Sal 


Ah? ~ - 
: Po «3 


(51) 
Equations (50) and (51) can be used jointly in a relaxation process, the 
only essentially new point being that when we estimate A from the quotient 
(51) the terms of order A‘ in the numerator would have to be computed from 
the A-value determined by the previous cycle. If h is small enough these 
terms will be small; they constitute the difference between the eigenvalues 
of (46) and (50) and we can easily show that this is second order with respect 
toh. We shall consider only the case ( — 0 at the end-points. Let A be 
the approximation to a given A of (44) found from (46) and let A® be that 
found from (50), nodal ¢ and q-values being similarly distinguished. If we 
multiply (46) by and (50) by G), subtract and sum over all nodal points 


we obtain (1) 9) 
> Id 


Vie | bh2{qi? }? Cy” ris) 0 


i > fp f212yr(by ree) 
and hence A) AW yh do \So ‘So - (52) 


172) 


y a Po oo So 


Further, for small enough /, (5) ~ ¢ for every nodal point and py — 0 by 
hypothesis so that ultimately A —- A™ for every mode and the result is 
true independently of / for the first (or one-signed) mode. 

This result is interesting in view of Bolton and Scoins’ result (48). If 
this latter were true under more general conditions it would mean that 
\ always corrects (in the right direction, since A@ is an underestimate of 
the true A. We can show that this is so and, moreover, that A® is always 
an improvement. 


7. Suppose that the range of integration, from 2 — 0 to 2 — 1 say, is 
divided into NV equal intervals each equal to h — 1/N. We shall consider 
only the case ¢(0) t(1) 0 and suppose q(7) to be bounded in (0, 1) 
and non-singular at the end-points. It is convenient to write (15) in 
terms of successive derivatives of ¢ at the point 0 rather than in terms of 
differences. If we then substitute in (44) we obtain the finite-difference 
relation 


¥ ro, ‘6 Y 2riv ¥ n¢ 
Apy So A "(Git ly 20) +90 S04 yh CO ! ghol*Cy" Trees (53) 
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holding at all internal points, from which we obtain the Rayleigh quotient 
by multiplying by (, and summing for all internal points, i.e. 

r > { NC, T ly 20) t % Cot jh? t shoh*Cs' 1 So 

> Po u 
If we suppose the higher derivatives to be expressed in terms of pivotal f- 
values and assume /& to be small enough for (53) to converge, (54) defines 


(54) 


Aas a function of the N—1 internal (-values, Its stationary values for all 
variations of the pivotal values, obtained by putting ¢A/ef, —— 0 for each f, 
give rise to (53) and are the true eigenvalues of (44). Neglecting terms in 
h? and higher powers of # in the numerator gives the quotient (47) and leads 
to the approximation A —— A“, In its full form the quotient (54) is equiva- 
lent to the extended form of Bolton and Scoins’ relation (48)+ and we must 
now examine the sign of the terms in A?. In the interval (say) 2 — 2, to 
x aw, we have, using trapezoidal summation, 

( C(ati(x) dx — bh(Sy Sy 4-0, Oy") + Of). (55) 
Summing (55) over each interval and remembering that { — 0 at the end- 
points we obtain, approximately, 


I 
WS Loh | Layer) dar [Lay — Saye") | + | (0a)? de. 


0 0 (56) 
The first term on the right hand of (56) vanishes by virtue of the boundary 
conditions since, by hypothesis, neither p(2) nor o(2) is singular at an end- 
point. We can therefore replace the 4?-sum in the numerator of (54) by 
the term , , 
ih | {C"(x)}2 dx = ph | Pe? dx, (57) 
0 0 
using the differential equation (44). Since p(x) > 0, by hypothesis, it 
follows at once that (48) is true for a variable p and that vy < 0, Moreover, 
replacing the integral (57) by its trapezoidal sum we can now write the 
quotient (54) as i 
>} NC, 1 Cy 20) Op lo t jyhPgi Co t O(h*) {Lo 
j > po Si 
and applying the variational condition to (58) leads to 
Cit ly (2—hPgot yh *qity = Oh), 


i.e. (50) is correct to a higher order of accuracy than (46). This result seems 


° 


A (58) 


+ Summing over all pivotal values effectively reduces each power of h in the numerator 
of (54) by one, but the sum in the denominator is of order 1/h, #0 the equivalence follows. 
Hoe 52 K k 
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particularly useful since if we first compute an eigenvalue and mode of 
(46) we can at once get an improved eigenvalue by substituting our com- 
puted mode in (51), i.e. without any further differencing of the solution, 
The mode can then be improved using (50) but subsequent changes are 
unlikely to alter the eigenvalue by very much. Incidentally we may notice 
in passing that although we have approximated to (49) by (50), this is not 
necessary. We can use (49) in conjunction with its proper Rayleigh 
quotient, which only involves replacing the factor 2--A‘gj, 12 in (51) by 


2 cos UN h h*qy. 


» 


Kxrample ~ 
To test these results we have applied them to an example previously 
solved by Fox (op. cit. 184) who uses it to illustrate the ‘difference-correc- 
tion’ method of standard finite differences. The problem is that of finding 

the mode of the Mathieu equation 
O° 4-(A—2cos 2r)t = 0 (59) 


which vanishes at x Qanda 7 and has no internal zero in the interval 


(0, a); this mode is symmetrical about 2 — 47. Fox uses an interval 


h iv and obtains a first approximation according to (46), i.e. in this case 


C,-+-C,3—(2—Ah?+- 2h? cos 225) Cy 0, (60) 
estimating A from the corresponding Rayleigh quotient, His final result for 
the eigenvalue is A? 00210 which, to this number of figures, is the 
correct eigenvalue of (60), This is then improved by differencing the corre 


sponding approximation for € and adding difference corrections according 


to (15), leading eventually to Ah? O-OL7T008, i.e. A 0-11029, which 
is known to be the correct eigenvalue of (59) to four decimals. The under- 
estimate given by the eigenvalue of (60) is clearly substantial in this case 
(almost 24 per cent, of the true value) and we have found that this result 
can be improved considerably if we use (50) instead. The eigenvalue 
obtained from this set of equations is Ah? 0-O1712, i.e. A O-1110 
and the corresponding eigen ector is given in Fig. 2. This lies between the 


eigenvector obtained from (60) and Fox's final corrected eigenvector and 
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although it is not nearly as accurate as the latter it leads, as we expect 
from the properties of the Rayleigh quotient, to a much improved eigen- 
value (now not much more than 4 per cent. in error). 


8. We can apply Allen and Southwell’s treatment of (1) to (45) except 
that there is not now a unique way of separating (45) into its two component 
equations. The simplest method is to write this equation in the form 

—- + 4q¢ A(a,y), vit. + 4qe A(x, y), (61) 
ox* cy 
where again gq — Ap—o. We can now approximate to the first equation of 
(61) along the x-line 301 by the equation 


c 
1 ¥ 

ne 1 890o Ay 

ox? 


which leads to the approximation 


4/ 
Ly 4 Ly — 2Lyeos| (Aqo)*h fat cos|(Aqo)!h]} — 0. (62) 
0 


Dealing similarly with the second of (61) along the y-line 402 we obtain 


4A, 


ly + C4 2f, cos] (Aqy)*h| 4 {1 —cos|(dqo)'A]} — 9, (63) 
) 


so that if we add (62) to (63) we get the approximation 
Cyt Ca 4 Cy 4+-Cy— 40, cos] (dq0)'h| 0. (64) 


This result can be closely represented, for small enough h, by the approxi- 
mation i 4-Setle+0q4—(4—h2 qo + hh*g?lo = 0 (65) 
and, as for (44), it differs by a term in h‘ from the standard finite-difference 
formula C4 et ls +%4—(4—A% Golly = 0. (66) 


Two other obvious ways of separating (45) into component equations are 
the arrangements 


> — Bax,y), “24+9t = —Biz,y) (67) 


cr 


9 2 
o* on 


s+ = Clz.y, —< C(x, 9). (68) 
cr cy” 


and 


Each of these gives a nodal approximation which is different from (64) and 
neither of these approximations is symmetrical with respect to nodal 
t-values, although we ought really to expect this symmetry on account of 
the symmetry of (45) with respect to the second derivatives. However, we 








500 S. C. R. DENNIS 


can combine the approximations obtained from (67) and (68) into the single 
approximation 
¥ i ‘ . > y . | ae a ] ; = > 
C4 Co 4 Cg + 4—40of2+ (h2qo—2)eo0s gi h}/{2+h¥qg—2cosqih} = 0 (69) 
and although (69) is not the same (or as simple) as (64) it leads to (65) if 
we expand in powers of / and retain only terms in h* and below. This lends 
more confidence to the use of (65) as a satisfactory approximation. 
The Rayleigh quotient corresponding to (65) is 
° e v ¥ 2 2 ¥9 
= | (Sit ty S31 So (4+ digi) loi So h oy CI 


- vo 
Zz Po So 


Ah? (70) 
which differs only by the term in #4 from the quotient corresponding to 
(66). If € — 0 on the boundaries of any rectangular region of integration 
we can show, using a method similar to that used for (44), that the corre- 
sponding form of (52) is now 
vo gan ES PPLE see 
z po sn Sy” 
where A”) now refers to an eigenvalue of (66) and A® to the corresponding 
eigenvalue of (65). Moreover, for the case of rectangular boundaries, 
Bolton and Scoins (5) have shown that (48) (A here being A” in our notation) 
is equally applicable to the two-dimensional Schrédinger equation and 
that v < Oif ¢ 0 on the boundaries. This is again interesting since (71) 
shows that, for small enough h, A > A” by a term of order h?. An in- 
vestigation along the lines of the previous section now proves much more 
difficult, however, and we shall not attempt it. Instead we shall suggest 
a somewhat tentative justification of the superiority of (65) by considering 
the results of two numerical examples. Both relate to the transverse vibra- 
tions of membranes, in which o = 0 in (45) and p represents the surface 
density of the membrane, ¢ being the transverse deflexion. 
Example 3 
The first relates to a rectangular membrane with sides in the ratio 2:3. 
The surface density is constant and the membrane is fixed (¢ = 0) along 
its edges. This case has previously been considered by Allen (6) as an 
example to illustrate the use of standard relaxation methods. We shall 
suppose the units to be non-dimensional, the length of the shorter side 
being unity. Taking the origin at one corner of the rectangle with the 
y-axis along the shorter side, the appropriate governing equation is 
V2Z+AZ = O (72 
and the exact solutions for the modes are 


oe 
c sin“ Jsin(nzy) (73) 
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with eigenvalues A = 7°({m?+-n?), (74) 
where m and n are any integers. The first mode (m = n = 1) possesses 
symmetry about both x = } and y = } and is the one which Allen has 
computed approximately. Using an interval h = 4 and the standard 


formula (66) (with g, =A) he obtains the final value A = 14-14, which 
underestimates the correct answer by rather less than 1 per cent. We 
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solved the same problem using (65) and our final answer on the same interval 
is A = 14-24; this is again an underestimate but of an amount little more 
than 0-1 per cent. The corresponding mode is given in Fig. 3 (the maximum 
detlexion being fixed at 1,000) and differs nowhere by more than a unit 
from the appropriate mode of (73) (when scaled to the same maximum 
deflexion). 
Example 4 
Our second example relates to a square membrane with a variable surface 
density, previously treated by Allen, Fox, Motz, and Southwell (7) by 
standard methods. With coordinate axes as before the appropriate non- 
dimensional governing equation is 
V24+ 2A(1—a—y+2ry)l = 0, (75) 


¢ vanishing on the unit square bounded by the axes and the lines x = 1, 
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y — 1. Using an interval h J;, Southwell et al. have computed the mode 
with no internal nodal lines using (66); their final accepted eigenvalue is 
A = 19-5661+ and they state that they would expect (from consideration 
of similar computations when the density is constant and the exact eigen- 
value known) that the true A should lie between 19-54 and 19-60. As our 
first test we substituted their final computed mode into (70) and found that 
it resulted in an estimate A — 19-62881. Residuals were then calculated 
from (65) and it was found necessary to modify the initial assumption for 
the mode by not more than one unit at a few nodal points to arrive at an 
acceptable solution, Substitution of this in (70) gave A — 19-62878, which 
is our final value. To see whether our value is an improvement in this case 
we have obtained an independent solution to the problem. 


9. The process we use is essentially that of the method of variation of 
parameters. The modes of (75) satisfy ¢ — 0 on the unit square and can 
therefore be expressed as an absolutely and uniformly convergent series 

f Y 3% a,sinizersin jry, (76) 
U 1 j l 


where ¢ and / are positive integers. If we substitute in the expression 


11 11 
U } | CV°C dady+-2A | | (l—ax—y- 2xry)l? dady (77) 

0 0 00 
and put ¢l ¢a,,, 0 for all integer values of m and n, i.e. we use the 


property that (75) is satisfied at a stationary value of (77), we obtain the 


doubly-intinite set of algebraic equations 


, , > ) 
}ar=(m*-+- n=) Ala 


64Amn ~ ~ f1—(—1)**Yf1—(— 1) + hay 
a fins dont (i? m*)*( )? n*)? 
(« #Am,j #n), (78) 

any equation of the set being obtained by assigning a pair of integer values 
to mand n. The characteristic solutions of (78) are the true modes of (75) 
expressed in the form (76); we cannot solve (78) exactly but we can approxi- 
mate to a given mode to any desired accuracy by standard relaxation 
methods. For example, we are looking for the mode with no internal nodal 
lines and a good approximation ignores all but the first term of (76), giving 
from (78) - 


(27?—A)a,, 0 


i.e. A 27°, with the corresponding mode 
r4 Sin max sin wy. 


Tt Subsequently amended to 19-564 by Miss G. Vaisey (see Southwell (8)); we worked in 
relation to the original paper cited. 
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We can improve this by taking account of the next important term in (76) 
(judged to be the one whose coefficient is associated with the next smallest 
factor on the left-hand side of (78)). In fact this is a. since it is obvious 
that (78) consists of two independent sets of equations involving firstly only 
those a,, for which i+-j is even and secondly only those for which i-+j is 
odd. Working to five decimals we then obtain the equations 

(27°—A)a,,—0-12978Aa,, = 0, 

0-12978Aa,,— (S27? —A)dgy = 0, 
from which we obtain a new estimate A = 19-62981 with the mode 

C = sinawrsin zy+0-0429 sin 27x sin 27y. 

By extending this process to sufficient equations we have computed the 
mode and eigenvalue to high accuracy, i.e. we have included all 
a,; ~ 0-00005 and taken into account enough equations to ensure that 
no a;,; we have omitted has any effect on computed values. The equations 
were solved by relaxation methods and the eigenvalue estimated from the 
corresponding Rayleigh quotient. It is easily shown that a,;; — a;,; and 


our values are given, correct to four decimals, in Table 2. 


TABLE 2. Values of a;; 





: I 2 3 4 s 6 
} 
I 1°0000 00016 o°0001 
2 0°0431 ooo12 0°0002 
3 o0016 00004 o°000!1 
4 o°0012 o*0001 
5 0°0001 o"o00!1 
6 0°0002 


The eigenvalue is A = 19-62884 and is the true eigenvalue of (75) to all 
the figures quoted. Returning to the relaxation solutions it follows that 
Southwell et al. are incorrect in their prediction of the bounds for A and that 
our value underestimates the true value by about three parts in a million, 
a result of quite surprising accuracy. We could add the general comment 
that, as far as we can say from these examples, there is a definite advantage 
in using Allen and Southwell’s type of approximation to solve equations of 
type (45). For the final part of our paper we shall return to equations with 
the first derivative present and make an extension of our previous results 
which is of more especial interest to ordinary differential equations, for 


reasons that we now explain. 


10. An extension of Numerov’s method 
The question of the desirability of using improved finite-difference 
formulae (i.e. with error terms of a higher order with respect to h) in the 
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solution of partial differential equations such as (1) or (45) is debatable. 
Such formulae require the introduction of more nodal points than are 
present in Fig. 1 and the relaxation process is more complicated. But for 
second-order ordinary equations in the normal form 


C"+-q(x)l = r(x) (79) 
a more accurate three-point formula is obtainable. It follows from (15) that 
A*(1 + 120°) Lo (0?+ 240°" a 31299 ++++)Eo (80) 


and we can use (79) to express 622) in terms of pivotal values of r—q{ so 
that a finite-difference equivalent of (79) is 

(1+ tyh?qy)0, + (1+ h?qg)S3— (2 — 8A? q0)So— tsh*(7y + 10r9 +15) 0 (81) 
with an error of only 5\,5°¢y. Equation (81) is generally ascribed to Numerov 
(9) but it cannot be applied in this way when the first-derivative term 
is present (as in (8)) because 6?f) would now involve 8?{, and more pivotal 
points must be introduced. Our previous results suggest how to overcome 
this. We first transform (8) to (21), apply Numerov’s method to this latter 
equation and then return to our original variable as before. This gives the 
pivotal approximation 

(1+ gh?) )C, + (1+ yh?bg) lye" —(2— Bh?dby)Sy— 
~phP(r, e+ L0rg +73 e") 0 (82) 

as an accurate three-point representation of (8). By the term ‘three-point’ 
here, of course, we mean involving only three pivotal values of ¢. If the 
functions ¢ and u are expressed analytically the formula is three-point in 
every sense; but if these functions must be expressed approximately in 
terms of pivotal values we shall need extra pivotal points (certainly for 
finding 4, and ¢,). This can always be done ab initio and does not alter the 
relaxation process at all. 


Exam ple » 


To test this result we have again used a problem solved by Fox (3, 73), 
i.e. to find an approximation to the solution of 


yn 


g + 


Dey 


‘+490+-1 0 (83) 
with end conditions ((0) = {(1) = 0. Using .n interval h }, Fox first 
obtains an approximation according to the standard formula 

(1+-$h)C, + (1— $h)f,—(2—4h225)fp +h? = 0 (84) 
which he then improves by calculating difference corrections and adding 
them in according to the full formulae (14) and (15). We give both his 
solutions in Table 3, the corrected final solution being not more than a unit 
out in the last decimal from the true solution. We have obtained this 
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accuracy directly using (82). We also give, for interest, solutions according 
to (12) and (23), although we do not expect either of these to give results 
comparable with (82). Allen and Southwell’s formula gives a better result 
than (23) but probably only, because the term involving p’ in ¢ (which 
the former fails to take account of) is identically zero in this case. In any 
case, both give better results than (84). 


TABLE 3 














¢(L. Fox) 4 C r4 
x Eqn. (84) Corrected Eqn. (82) Eqn. (12) Eqn. (23) 
oo o"o oa oo o'o oo 
0°25 01236 o1218 o1218 O1213 o°122. 
"50 O1S573 o1s4! O1s4i or1§40 o°1557 
O75 O1104 O1075 0'1076 0" 1080 O"1094 
1°00 ° ° ° ° ° 











11. Conclusions 

The investigation of this paper was undertaken to discover why Allen 
and Southwell’s method of approximating to (1) should be superior to stan- 
dard approximations. There can be little doubt that it is, for it deals more 
effectively with the terms involving first derivatives. But we have shown 
that (1) is rather a special case since the coefficients of the first-derivative 
terms are related so that a term which may otherwise be important cancels 
identically. That this term can be important is clearly indicated in example 
1; the investigation here relates to the ordinary differential equation (8) 
but the essential features are the same. In fact, cancellation of the term 
in question (arising from p’) is not possible for ordinary differential equa- 
tions unless as is exceptional, p(2) is constant, as inexample 5. The investi- 
gation has led us to a number of alternative approximations which, as far 
as we can see, possess all the features of the removal of the first-derivative 
term from the equation, i.e. they behave with respect to convergence as if 
this term were absent (a case in which it has always been known that finite- 
difference methods give more satisfactory results). In the case of ordinary 
differential equations we can extend the basic method of approximation 
to take account of higher differences without introducing more pivotal 
points (an extension previously limited to equations with the first-derivative 
term absent). 

The central section of the paper, which deals with the eigenvalue prob- 
lems governed by (44) and (45), involves questions of a rather different 
nature since now no first-derivative terms are present in the equations. 
The essential point here is whether Allen and Southwell’s process still 
gives improved results. Subject to certain limitations (some of which could, 
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perhaps, be removed), we have shown that the approximation to (44) is 
certainly improved, The fact that the error term in (50) is of order h® is 
of some interest in view of the fact that the derivation of the approximation 
(12), from which it is ultimately obtained, does not require ¢ to possess 
derivatives beyond the second order. This is distinct from the improved 
order of accuracy which would result from applying Numerov’s approxima- 
tion (81) since this requires higher derivatives of ¢ to exist, although the 
distinction is a theoretical rather than a practical one. In the case of the 
partial differential equation (45) our investigation is more tentative. Even 
if we restrict ourselves to the case € — 0 on a rectangular boundary, it is 
more difficult to prove that (65) is a more accurate approximation than (66). 
The evidence of the numerical example suggests, however, that it is. 
Consider, for example, the problem relating to the transverse deflexion of a 
uniform rectangular membrane (o¢ 0, p 1 in (45)). In this case (71) 


reduces at once to A(2 AD + LhAf A(2}2, 
which we can write in the form 
\? AD 4 BAAD 2 ORS). (85) 


In the particular examplet (cf. Southwell (8), 383) of a uniform mem- 
brane bounded by a unit square it is known that, for the fundamental mode, 


\" 19-72347 when 4 — 3. Neglecting terms in h* in (85) we obtain 
VG 19-73930 as against the exact value A = 27? 19-73921, to five 


decimals. Again, if we consider the membrane with variable density p(2, y) 
and suppose that in (71) we can approximate to the pivotal p-values by 
some mean value f (since, in effect, the right-hand side of (71) isan average), 
then we can approximate to A® by the relation 

A® = AM+4 LAA5{AM}3, (86) 


An obvious choice for 6 would be the integrated mean value, i.e. 


| ( pla, y) dady 
Db 

| | dxdy 

D 
where PD is the domain bounded by the membrane. In example 4, for 
instance, this gives p 1 and, using Miss Vaisey’s value A = 19-564 
when h fy. we obtain A® — 19-626 from (86), a very great improvement 
obtained with no substantial additional labour. 

In general the use of Allen and Southwell’s type of approximation raises 
very wide issues, of course, e.g. when we may best expect a function to be 


+ IT am indebted to a referee for pointing out this example. 
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represented more adequately by expansions other than polynomials. Such 


general questions cannot readily be answered and this paper provides only 


a limited amount of evidence that may be worth consideration in specific 


cases, 
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SUMMARY 
This paper describes a method of solving an algebraic equation by determining 
high powers of an associated matrix using the Cayley- Hamilton theorem. The method 
described will be found to have practical applications to linear equation analogue 
computers for finding the zeros of polynomials. 


1. Introduction 

ALTHOUGH various numerical methods are available in literature (1, 2, 3) 
for finding the zeros of polynomials. not much attention has been paid so 
far to the application of matrix meu ods of solution (viz. converting the 
problem to an eigenvalue problem) which has been suggested by Aitken 
(4). Also the methods of finding zeros of polynomials in all the analogue 
computers constructed so far have been direct (5), using non-linear and 
trigonometric function potentiometers. It is the object of this paper to 
point out that Aitken’s method (4) can be readily mechanized with great 
advantage either in linear equation computers in which matrix multiplica- 
tions can be accomplished (Barker (6)), or in direct secular-equation 
solvers (5). This, together with the acceleration process herein suggested, 
results in the significant improvement of the method of determination 
of real or complex roots of equations, with either real or complex 
coefficients. 


2. The matrix method 
It is well known (7) that any polynomial of degree n 


ay+a,xr+a,z?+...+a,_,27" "+2" = 0 (1) 


2 


has the associated companion matrix (n x 7) 


A oe . 2.4% 8 Ay 
lo... 0 —a, (2) 
oe OS 4S le -An-1- 


Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 4, 1960) 

















aa ae 
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whose eigenvalues are the roots of (1). Hence, the problem reduces to 
the determination of eigenvalues of A for which a variety of methods are 
already available (1, 2, 3, 8). Among these, the matrix power method 
is the one that is the easiest and quickest, for the determination of the 
eigenvalue of A with largest modulus. 

The matrix power method essentially consists in repeatedly multiplying 
an arbitrary column vector X by the matrix A so that the recurrence 
relation AXom-) — Xim) (3) 
is satisfied. 

For purposes of computation it is desirable to reduce the largest of the 
components of X’” to unity and each time this normalized vector is fed 
back for getting the new iterated vector X’*), Therefore when this 
process converges (to a desired accuracy) at the end of m cycles, we have 


AX™ — A Xim+d, (4) 


where X is almost identical with X°"*” and A, is the root with largest 
modulus. This method, as is easily seen, thus involves the successive 
multiplications by the matrix A, say m times. Even this labour could be 
considerably reduced by a technique suggested in the next section. 

Once the largest eigenvalue is known, the other eigenvalues are obtained 
either by using, (1) the deflation method, (2) Richardson’s purification 
process, (3) the Gram—Schmidt process described in the references cited 
above, or by reducing the polynomial itself to one of a lower degree by 
synthetic division and again solving the ‘reduced equation’ by applying 
the matrix power method. However, the rounding-off errors in these 
processes have to be carefully considered when high accuracy is aimed at. 

Once all the roots are obtained to a sufficient approximation, they can 
be further refined either by the Bairstow—Lin process (3) or Jahn’s (9) 
method: 


3. The method of acceleration 


Since we know by the Cayley-Hamilton theorem that any matrix A 
should satisfy its own characteristic equation, we have by (1) 


—(a,1+a,A+...+a,_,A"-!) = A". (5) 


Therefore any higher power of A, say A” (m > n) is always expressible 
in terms of a matrix polynomial in A, involving only the first (n—1) 
powers of A. So by choosing m to be large, quick convergence can be 
obtained to the eigenvector X, corresponding to the largest of the eigen- 
values A,, just by a single multiplication of an arbitrary vector by A”, 
which is known, provided an initial algebraic calculation is made to 
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express A” in terms of A’ [where r = 1, 2...., (n—1)]. When the eigen- 
vector X, is thus known, a single multiplication by A could yield A, also. 

The computation of A’ does not present serious difficulties because of 
the inherent nature of the companion matrix A (even though it may be 
difficult for any general matrix). The peculiarity with A is that all the 
first (n—1) columns of A? are the same as the second to the nth columns 
of A, all the first (n—1) columns of A® are the same as the second to the 
nth columns of A?, and so on. Consequently it is necessary to determine 
only (n—2) columns for getting A’ (r = 2,3,...,n—1). Thus, even if no 
analogue computer is available, these multiplications could be carried out 
with great ease, using desk calculators. 


4. Example 
These ideas will be clearer if a simple example is considered, namely, 
the determination of the largest root of the familiar equation 


a—z—] 0. (6) 
Here A 0 0 1), A2 0 1 07, 
lL Oo ] oO 1 i 
0 1 0 ; oe 1 


and A? — A--I by (5). Therefore, 
Als — [(A3)8]? = 37A2+ 49A-+ 281 

28 37 49 

49 65 86 

37 49 65 
Taking a general arbitrary vector X® = (1,1, 1), 

ABX(0 200 (0°570, 1-000, 0-755) 200 (X“)), 

Therefore AX8 1-325 (0-570, 1-000, 0-755). 


The value of A, is 1-325 (the exact value accurate to 4 figures . . . 1-3247). 
If still higher powers of A are used much higher accuracy can be obtained. 


5. Conclusion 

The method described above has several limitations. If the eigenvalues 
of (2) are all distinct and real, there is no difficulty in applying the method. 
However, it is possible that the dominant eigenvalues are either (a) nearly 
equal, or (b) complex pairs (even if A is real), or (¢) roots of equal modulus 
(real or complex), or (d) repeated (real or complex) roots. 

For the case (a), the convergence can be obtained without any difficulty 
taking a high enough power of A. 
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But in the case of (+), the matrix multiplication of an arbitrary real 
vector X would not lead to a convergent sequence (even if each pair is 
well separated from the next root (real or complex) and the elements of 
iterates would oscillate instead of tending to a limit. From a knowledge 
of three successive iterates (in the case of well separated roots) or equally 
spaced iterates (if several roots are quite close), e.g. A™X®, A™-X®, 
and A”**4X where A is any properly chosen integer, it is possible to 
determine the dominant complex pair of roots (10). Instead of having 
an arbitrary real vector, if a complex arbitrary vector is chosen, the pair 
of roots can be obtained directly. In such a case, however, the largest 
component of iterated vector must each time be normalized to (1+-7.0) 
and iterations performed in the usual way. 

For the case (c) when there are real roots of equal modulus the situation 
is practically analogous to case (b), but the successive iterates AX, 
A”™!X®, ete., will be alternately parallel to one of the two eigenvectors 
and hence the multiplicity will be revealed. In this case it is always 
possible to add a constant real matrix pI to the matrix A, thus altering 
the moduli of the real roots from -+-A to (+A+), and hence separating 
them. The same method can also be applied to separate real and complex 
roots having the same moduli. 

For the case (d), there are certain difficulties of convergence which are 
intimately connected with the properties of matrix (2). If the matrix (2) 
is diagonalizable, the repeated roots will have an independent set of 
eigenvectors associated with each one of them, and taking any general 
arbitrary vector the convergence will not be on one of the eigenvectors, 
but will only be a vector which is a linear combination of both these 
vectors. However, there will not be any difficulty of convergence in such 
a case, 

But in the case of repeated eigenvalues (real or complex), it can happen 
that the number of linearly independent vectors associated with this 
eigenvalue is less than the multiplicity. In this case the matrix cannot 
be diagonalized and the convergence to a limit will be extremely slow. 
Unfortunately there are no simple and direct methods of testing the 
diagonalizability of a matrix. In such a case, therefore, one has either to 
go to very high powers of the matrix A and use Aitken’s 5?-process (3) 
or to use a method developed by the author which is akin to Graeffe’s 
technique. The latter method will be published elsewhere. 

The one great advantage of the matrix method is that it breaks up the 
polynomial into a set of linear equations, so that a linear equation analogue 
computer would serve the purpose better than a polynomial com- 
puter, which involves a number of non-linear computing elements (5). 
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